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ABSTRACT 

Radiation feedback is typically implemented using subgrid recipes in hydrodynamical 
simulations of galaxies. Very little work has so far been performed using radiation- 
hydrodynamics (RHD), and there is no consensus on the importance of radiation 
feedback in galaxy evolution. We present RHD simulations of isolated galaxy disks of 
different masses with a resolution of 18 pc. Besides accounting for supernova feedback, 
our simulations are the first galaxy-scale simulations to include RHD treatments of 
photo-ionisation heating and radiation pressure, from both direct optical/UV radia¬ 
tion and multi-scattered, re-processed infrared (IR) radiation. Photo-heating smooths 
and thickens the disks and suppresses star formation about as much as the inclusion of 
(“thermal dump”) supernova feedback does. These effects decrease with galaxy mass 
and are mainly due to the prevention of the formation of dense clouds, as opposed to 
their destruction. Radiation pressure, whether from direct or IR radiation, has little 
effect, but for the IR radiation we show that its impact is limited by our inability to 
resolve the high optical depths for which multi-scattering becomes important. While 
artificially boosting the IR optical depths does reduce the star formation, it does so 
by smoothing the gas rather than by generating stronger outflows. We conclude that 
although higher-resolution simulations, and potentially also different supernova imple¬ 
mentations, are needed for confirmation, our findings suggest that radiation feedback 
is more gentle and less effective than is often assumed in subgrid prescriptions. 
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1 INTRODUCTION 

To first order, gravity describes the formation of structure in 
the Universe ( Peebles V Yu|1970l|Zerdovich|1970 ). The for¬ 
mation of galaxies in dark matter (DM) halos also requires 
radiative cooling to relieve pressure and dissipate angular 
momentum ( Binney|1977 Rees V Ostriker|1977 Silk|1977 ). 
It is also well established that in order to halt the collapse of 
gas into galaxies, dense substructures, and eventually stars. 


counteracting feedback processes are required (e.g. White V 
Rees 1978). Without feedback, galaxies collapse and form 


stars too efficiently, compared to observations. 

Early simulations focused on feedback in the form of 
supernovae (SNe; e.g. Katz I992| Navarro V White 11993 ) 


and later active galactic nuclei (e.g. Di Matteo et al.||20Q5 
Booth V Schaye|20Q9 Dubois et al.|20I(5 ), where the latter 

is thought to be dominant in massive (“L > L*”) galax¬ 
ies (Bower et al. 2006). However, simulations that include 


those feedback processes still struggle to produce galaxies 
that match observations in terms of their star formation 
histories and morphology, ( Scannapieco et al.|20I2 ). 


Analytical work by e.g. Thompson et al. (2005) and 


Murray et al. (2005 2010 2011) suggests that radiation 


feedback may be an important missing ingredient. Recent 
hydrodynamical simulations therefore often enlist stellar ra- 


diation in their subgrid feedback models (e.g. 

Brook et al. 

2012 

Agertz et al.||20I3| Stinson et al.||20I3| | 

loskar et al. 

2014 

Ceverino et al.||20I41 |Kannan et al.||20I4a|bl [Hopkins 

et al. 

20I4| Agertz & Kravtsov||20I5|). The added radiation 


E-mail: joki@strw.leidenuniv.nl 


feedback usually contributes directly to direct suppression 
of star formation, and increases galactic outflows, which 
can expel the gas altogether and enrich the intergalactic 
medium (IGM) with metals. The idea of radiation feedback 
has proven so successful that most cosmological simulations 
nowadays invoke it in some form, although the implemen¬ 
tations vary a lot, and they are often motivated empirically 
rather than physically. Radiation feedback on galactic scales 
is usually modelled with subgrid recipes in otherwise purely 
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hydro dynamical (HD) codes. These HD recipes must make 
a number of assumptions about e.g. the absorption of pho¬ 
tons, mean free paths, and shielding. They can thus only to 
a limited degree be used to investigate how important radi¬ 
ation is for the formation and evolution of galaxies, and how 
the radiation interacts with the baryons, i.e. how radiation 
feedback actually works. 

The recent literature on simulations of galaxy evolu¬ 
tion usually considers three radiation feedback processes: 
photoionisation heating of gas, direct pressure from ionis¬ 
ing photons, and indirect pressure from reprocessed, multi¬ 
scattering, infrared (IR) photons. Simulations often contain 
only a subset of these processes, and there is no general con¬ 
sensus on the importance of radiation feedback as a whole, 
or on which of these processes dominate under which cir¬ 
cumstances (see {4.4). 

A more assumption-free and physically correct descrip¬ 
tion of radiation feedback requires the use of radiation- 
hydrodynamics (RHD), which models the emission and 
propagation of photons and their interaction with the gas 
self-consistently. RHD can help tell us if and how radiation 
feedback works, and this information can then be used to 
improve HD subgrid recipes of radiation feedback. 

However, RHD is both complex and costly compared 
to HD. For the most part, it has therefore not been used 
directly in simulations of structure formation, or more gen¬ 
erally in studies of galaxy evolution. In recent years however, 
the use of RHD has been on the rise in computational astron¬ 
omy, and RHD implementations have evolved towards being 
usable in cosmological and galaxy-scale simulations that re¬ 
solve the interstellar medium (ISM) ( [Wise et al.||2012b|a| 
Pawlik et al.|2013 Wise et al.|2014 ) 


In Rosdahl et al. ( 2Q13[ hereafter R13) we presented 
an RHD implementation in the cosmological code RAMSES 
( Teyssier]|2002 ), which we called RAMSES-RT. In that paper 
we modelled the emission and propagation of photons and 
their interaction with hydrogen and helium via ionisation 
and heating. In Rosdahl & Teyssier (2015 hereafter R14), 
we added two aforementioned processes to the implementa¬ 
tion, which are though to be relevant for galactic feedback: 
radiation pressure, i.e. momentum transfer from photons to 
gas, and the diffusion and trapping of multi-scattered IR 
radiation in optically thick gas. 

In the present paper, we use the RHD implementation 
that we have detailed in the two previous papers to study the 
effect of stellar radiation feedback on galactic scales. We use 
a set of RAMSES-RT simulations of isolated galactic disk sim¬ 
ulations, where we include stellar radiation feedback, com¬ 
bined with “thermal dump” SN feedback. The main ques¬ 
tions we attempt to answer are: 


• What role does stellar radiation feedback play in reg¬ 
ulating galaxy evolution, and how does this role vary with 
the mass and metallicity of the galaxy? 

• How does the interplay of radiation and SN feedback 
work? Specifically, does radiation boost the effect of SNe? 

• Where stellar radiation feedback plays a role, what is 
the dominant physical process: photoionisation heating, di¬ 
rect pressure from the ionising photons on the gas, or in¬ 
direct pressure via dust particles UV and reprocessed IR 
radiation? 


In this paper, we study the effects of turning on the 


stellar radiation in galaxies, while making minimal assump¬ 
tions about what happens on unresolved scales. While using 
RHD implies radiation feedback is modelled from “first prin¬ 
ciples” , we stress that it is still necessary to make a number 
of approximations, both in the modelling of the radiation 
itself and in its interaction with gas and dust. Also, and 
importantly, although we resolve the ISM to some extent, 
we do not resolve molecular clouds, the scales at which the 
radiation feedback originates, and at which the radiation 
couples most efficiently with the gas. We expect the current 
simulations to give us hints as to what radiation feedback 
does in reality, and, equally importantly, to teach us what 
improvements in modelling and resolution are required in 
future work. 

The structure of the paper is as follows. In ^we present 
an overview of the code, the setup of galaxy disks of three 
masses, and details of the modelling of gas, stellar popu¬ 
lations, and feedback. In we present the results, where 
we successively incorporate SN and radiation feedback pro¬ 
cesses and compare their effects on the galaxies. We focus 
on the suppression of star formation and the generation of 
outflows, study how radiation feedback plays a role, and ex¬ 
amine trends with galaxy mass and metallicity. In Q we 
discuss and justify our main findings on analytic grounds, 
demonstrate how they are limited by resolution, probe what 
effects we can expect when the resolution is increased be¬ 
yond the current limits, and qualitatively compare our re¬ 
sults to previous publications. Finally, in ^we summarise 
our main conclusions and discuss interesting future direc¬ 
tions. The appendices provide details on the model we use 
for stellar population specific luminosities and convergence 
tests. 


2 SIMULATIONS 


We use RAMSES-RT (R13, R14), an RHD extension of the 
adaptive mesh refinement (AMR) code RAMSES (Teyssier 


2002). RAMSES models the interaction of dark matter, stel¬ 


lar populations and baryonic gas, via gravity, hydrodynam¬ 
ics and radiative cooling. The gas evolution is computed us¬ 
ing a second-order Godunov scheme for the Euler equations, 
while trajectories of collisionless DM and stellar particles are 
computed using a particle-mesh solver. RAMSES-RT adds the 
propagation of photons and their on-the-fly interaction with 
hydrogen and helium via photoionisation, heating, and mo¬ 
mentum transfer; and with dust particles via heating and 
momentum transfer. The code solves the advection of pho¬ 
tons between grid cells with a first order moment method 
and closes the set of radiation transport equations with the 
Ml relation for the Eddington tensor. The trapped/stream¬ 
ing photon scheme presented in R14 describes the diffusion 
of multi-scattering IR radiation. The radiation in a photon 
group, defined by a frequency interval, is described in each 
grid cell, by the radiation energy density E (energy per unit 
volume) and the bulk radiation flux F (energy per unit area 
per unit time), which corresponds approximately to the ra¬ 
diation intensity integrated over all solid angles. RAMSES-RT 
solves the non-equilibrium evolution of the ionisation frac¬ 
tions of hydrogen and helium, along with photon fluxes and 
the gas temperature in each grid cell. 

Because the timestep length, and therefore the com- 
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Table 1. Simulation parameters for the three disk galaxies. The listed parameters are, from left to right: Galaxy acronym used throughout 
the paper, r^circ^ NFW circular velocity, for the IC generation, i^vir^ halo virial radius (defined as the radius at which the DM density is 
200 times the critical density at redshift zero). I/box' simulation box length, Mhaio- DM halo mass, Mbisk- disk galaxy mass in baryons 
(stars+gas), fgas' disk gas fraction in the ICs, Mbuige- stellar bulge mass in the ICs, A^part: Number of DM/stellar particles in the ICs, 
m*: mass of stellar particles formed during the simulations, Axmax: coarsest cell resolution, finest cell resolution, ^disk^ disk 

metallicity. 


Galaxy 

acronym 

'Ccirc 

[km s“^] 

Rvir 

[kpc] 

hbox 

[kpc] 

-^halo 

[Mo] 

-Mdisk 

[Mq] 

/gas 

M^bulge 

[Mol 

Npart 

[Me] 

AXmax 

[kpc] 

^^min 

[pc] 

^disk 

[Z©] 

g8 

30 

41 

150 

IQio 

3.5 X 10® 

0.5 

3.5 X 10^ 

10® 

600 

2.3 

18 

0.1 

g9 

65 

89 

300 

10“ 

3.5 X 10® 

0.5 

3.5 X 10® 

10® 

600 

2.3 

18 

0.1 

GlO 

140 

192 

600 

10“ 

3.5 X lOio 

0.3 

3.5 X 10® 

10® 

10^ 

4.7 

36 

1 


putational load, scales inversely with the speed of light c, 
we apply the so-called reduced speed of light approximation 
( Gnedin Abel||2001[ R13) in runs that include radiation, 
to maintain a manageable computing time. In this work, we 
use a light speed fraction fc = 1/200, i.e. free-streaming 
photons are propagated at a speed c = c/200, such that the 
timestep is most of the time limited by non-RT conditions, 
and the slow-down due to RT is only about a factor 2-3 
compared to HD simulations, depending on the number of 
photon groups and processes included (and the inclusion of 
SN feedback, which limits the timestep as well). We showed 
in R13 that larger values for fc than we have chosen here are 
preferable in simulations of galaxy evolution in order to ac¬ 
curately capture the expansion speed of ionisation fronts in 
the ISM, but the light speed convergence tests presented in 
Appendixj^indicate that our results are robust with respect 
to the chosen light speed. 


We run simulations of isolated rotating disk galaxies 
of baryonic mass 3.5 x (10®, 10^, 10^°) Mq consisting of gas 
and stars embedded in DM halos of masses 10^°, 10^^, 
and 10^^ Mq, respectively. The simulation sets, named g8, 
g 9, and GlO, after the order-of-magnitude of the baryonic 
masses, are presented in Table and the parameters listed 
in the table are explained in what follows. The baryonic mass 
of the most massive galaxy (gIO) is comparable to that of 
the present-day Milky-Way (MW). 


2.1 Initial conditions 


The initial conditions (ICs) are gener ated with the 
Mak eDisi^ cod e by Volker Springel (see Springel et al. 
2005 Kim et ^|2014 ). The DM halos follow an NFW den¬ 
sity profile ( Navarro et al.|1997 ) with concentration param¬ 
eter c = 10 and spin parameter A = 0.04. We model the 
dark matter in each halo with A^part collisionless particles of 
identical mass. The initial disk consists of gas cells and A^part 
identical mass stellar particles, both set up with density pro¬ 
files that are exponential in radius and Gaussian in height 
above the mid-plane. The galaxies also contain stellar bulges 
with mass one tenth of the stellar disk mass, represented by 
0.1 A^part particles. The stellar particles that are present at 
the beginning of the simulation do not perform any feed¬ 
back. The initial gas profiles do not enforce exact hydrostatic 
equilibrium. However, the initial (few million years) stabili¬ 
sation of the galaxy, which manifests itself in contraction of 
the inner dense gas and expansion of the outer diffuse gas, 
is minor, as can be inferred from plots of the star formation 
rate (e.g. Fig.|^. The initial temperature of the gas disk is 
T = 10^ K, and the disk metallicity, ^disk, is set to a con¬ 
stant value, either 0.1 or 1 times Solar (see Table[^, with the 
metal mass fraction in the Sun taken to be Zq — 0.02. The 
circumgalactic medium (CGM) initially consists of a homo¬ 
geneous hot and diffuse gas, with nn = 10“® cm“®, T = 10® 
K and zero metallicity. The cutoffs for the disk’s radial and 
vertical gas profiles, which mark the transition between the 
disk and GGM, are chosen to minimize the density contrast 
between the disk edges and the GGM. 


For g8 and g 9, the host DM halos are disproportion- 
ally low in mass, compared to results from abundance¬ 
matching ( Moster et al.|20'T3 ) and cosmological simulations 


that match the observed galaxy mass function (Schaye et al. 


2015). These under-massive DM halos are not a major issue 


for the current work, however. We are primarily interested 
in comparing the relative effects of different feedback pro¬ 
cesses on the properties of the galaxy disk, for which the 
dark matter profile does not play an important role. To ver¬ 
ify that our results are insensitive to the mass of the host 
halo, we have run counterparts of the least massive galaxy, 
g8, with the halo mass increased to a more realistic value 
Mhaio = 7 X 10^® Mq (i.e. an increase by a factor of seven 
compared to Table [^, while keeping the same resolution. 
We confirmed that while the simulations were more expen¬ 
sive due to the increased size of the box and number of DM 
particles, the results were not affected. 


2.2 Star formation 


Star formation follows a standard Schmidt law. In each cell 
where the gas density exceeds the chosen star formation 
threshold 


n* = 10 cm“®. 

(1) 

gas is converted into stars at a rate 


— Cffp/tff, 

(2) 


where p is the gas density and Cff = 0.02 is the star for¬ 
mation efficiency per free fall time, tff = [37r/(32Gp)]^^^, 
where G is the gravitational constant. Gollisionless parti¬ 
cles of mass m*, representing stellar populations, are formed 


^ Adapted to generate RAMSES-readable format by Romain 
Teyssier and Damien Chapon. 
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stochastically from the gas, with the probability of forming 
one drawn from a Poissonian distribution (for details see 
). Table lists the stellar particle 

masses used in the simulations. In addition to the density 
threshold for star formation, we also do not allow stars to 
form in gas warmer than Tj— 3000 K, where /x is the 
average particle mass in units of the proton mass. We note, 
however, that our results are insensitive to increasing or even 
removing the temperature threshold. 


Rasera & Teyssier 2006 


2.3 Supernova feedback 


We model SN feedback with a single injection from each 
stellar particle into its host cell, 5 Myr after the parti¬ 
cle’s birth, of mass mej = ?7 sn x m*, and thermal energy 
esN = TjSN X 10^^ e rg m,./10 Mq. W e use ? 7 sn = 0.2, roughly 
corresponding to a Chabrier] ( |2003 ) stellar initial mass func¬ 
tion (IMF). We neglect the metal yield associated with stel¬ 
lar populations, i.e. the stellar particles inject zero metals 
into the gas. 

At our resolution, the “thermal dump” SN feedback 
model that we use is known to suffer from numerical over¬ 


cooling (e.g. Creasey et al. 2011 Dalla Vecchia & Schaye 
2012 Creasey et al.||2013 ), but we use it here, because it is 
simple and because it allows us to investigate how far radi¬ 
ation feedback can go to compensate for its low efficiency. 
The coupling between radiation and SN feedback, which we 
study in §3.1.4[ could depend on the choice of SN feedback 
model. More efficient SN feedback might either be amplified 
more efficiently by the stellar radiation to suppress star for¬ 
mation and increase outflow rates, or conversely, it might 
dominate completely over the effects of radiation feedback 
and render it negligible. These considerations are beyond 
the scope of the present paper, but in future work we will 
combine radiation feedback with more efficient recipes for 
SN feedback, to find what combination produces best agree¬ 
ment with observations (which is not the point of this paper) 
and to study how the interplay of the feedback processes is 
affected. 


2.4 Gas thermochemistry 


We evolve the thermochemistry semi-implicitly with the 
method presented in R13. The method tracks the non¬ 
equilibrium cooling rates of hydrogen and helium, here as¬ 
suming zero incoming photon flux. The ionisation fractions 
of hydrogen and helium are stored in each cell as three pas¬ 
sive scalars, which are advected with the gas. We assume hy¬ 
drogen and helium mass fractions X = 0.76 and Y = 0.24, 
respectively, and Solar ratios for the metal species, i.e. we 
track a single scalar representing the metal mass fraction in 
each cell. 

We add the contribution from metals to the cooling rate 


using tables generated with CLOUDY (Ferland et al. 1998), 


assuming photoionisation equilibrium with the redshift zero 


Haardt X Madau (1996) UV background. With metal cool¬ 
ing, the gas can in principle cool non-adiabatically to ~ 10 
K. We do not model the change in the metal cooling rate 
with the local radiation flux, which may affect galaxy evolu¬ 
tion (e.g. Cantalu^|2010 Kannan et al.|[^14b ). In future 
work, we will consider more realistic metal cooling, which 
takes the local radiation flux into account. 


2.5 Adaptive refinement 

In the adaptive refinement scheme of RAMSES, cells can be 
split into 8 child cells of width half that of the parent. The 
width of a cell is determined by its refinement hierarchy 
level by Axi = Lbox/ 2 ^, where Lbox is the simulation 
box width. The maximum and minimum cell widths, Axmax 
and Axmin, are determined by the enforced minimum and 
maximum allowed refinement levels in a simulation, which 
in this work are Axmax = 2 — 5 kpc and Axmin = 18 — 36 pc, 
depending on the simulation set (see Table [^. Adaptive re¬ 
finement follows mass: a cell is refined if it contains 8 or more 
collisionless particles, if the cell gas mass mceii > 12 m*, or 
if Ax is more than a quarter of the local Jeans length. 


2.6 Artificial “Jeans pressure” 


We impose a pressure floor on gas to prevent artificial frag¬ 
mentation below the Jeans scale (Truelove et al. ]T9^ . The 
Jeans length scale for a self-gravitating cloud is 



( 3 ) 


where Cg = ^y/ceT/mp is the sound speed, and we as¬ 
sumed a ratio of specific heats of 7 = 1.4, appropriate for 
a monatomic gas. From Eq. the requirement that the 
Jeans length is resolved by at least N cell widths becomes a 
temperature floor of the form 

— > ( NAx y 

1 K "" 1 cm“^ \16 pcy 

We apply this floor in the form of an effective temperature 
function. 


Tj = To nn/'R*, (5) 

where we use To = 500 K in all our simulations, ensuring 
that the Jeans length is resolved by a minimum number 
of 6 cell widths in g8 and g 9 and 3 cell widths in GlO. 
The pressure floor is non-thermal, and added to the physical 
temperature, T, and hence we can have T ^ Tj. 


2.7 Radiation feedback 


We include the emission and propagation of stellar radi¬ 
ation, and its interaction with the gas. The mass-, age- 
and metallicity-dependent stellar specific luminosities are 
extracted on the fly from the spectral energy distribution 


(SED) model of Bruzual & Chariot (2003), as described in 
R13, assuming a Chabrier (2003) IME. Stellar particles in¬ 


ject photons into their host grid cells at every fine RHD 
timestep. 

We bin the radiation into five photon groups, defined by 
the photon energy intervals listed in TableThe groups are, 
in order of increasing energy, IR, optical, and three groups 
of ionising ultraviolet (UV) photons, bracketed by the ion¬ 
isation energies for Hi, Hel, and Hell. We include the first 
two groups only in runs with radiation-dust interactions, 
while we include the three UV groups in all runs with ra¬ 
diation. Appendix [A| describes how the stellar luminosities, 
photon group energies, and ionisation cross sections are de¬ 
rived from the SED model. Table lists typical values for 
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Table 2. Photon group energy (frequency) intervals and properties. The energy intervals defined by the groups are indicated in units of 
eV by eo and ei (in units of Angstrom by Aq and Ai). The next four columns show photon properties derived every 5 coarse time-steps 
from the stellar luminosity weighted SED model (see Fig. |A1| and surrounding text). These properties evolve over time as the stellar 
populations age, and the approximate variation is indicated in the column headers, e denote the photon energies, while (Jho CTHeo cind 
(jHeii denote the cross sections for ionisation of hydrogen and helium, respectively, k is the dust opacity. The gas opacity scales with the 
gas metallicity, Ki = ki where i denotes the photon group. 


Photon 

group 

eo [eV] 

ei [eV] 

Ao [A] 

Al [A] 

e [eV] 
±10% 

(Thi [cm^] 

±5% 

(THei [cm^] 

±5% 

(THeii [cm^] 

±5% 

k [cm2 g ij 

IR 

0.10 

1.00 

1.2 X 10^ 

1.2 X 10^ 

0.6 

0 

0 

0 

10 

Opt 

1.00 

13.60 

1.2 X 10"^ 

9.1 X 10^ 

5.5 

0 

0 

0 

10® 

UVhi 

13.60 

24.59 

9.1 X 102 

5.0 X 102 

18.0 

3.3 X 10-1® 

0 

0 

10® 

UVHel 

24.59 

54.42 

5.0 X 10^ 

2.3 X 10^ 

33.4 

6.3 X 10-19 

4.8 X 10-1® 

0 

10® 

UVHell 

54.42 

oo 

2.3 X 102 

0 

60.0 

9.9 X 10-20 

1.4 X 10-1® 

1.3 X 10-1® 

10® 


the energies and cross sections, along with their variations 
over the simulation run-time. 

An important advantage of the moment method that 
we use for the radiative transfer is that the computational 
cost, i.e. the runtime of the simulations, is independent of 
the number of radiation sources. With the alternative class 
of ray-tracing methods (e.g. Wise Abel|20lT ), the compu¬ 
tational cost increases more or less linearly with the number 
of sources, which requires remedies to keep down the com¬ 


puting cost, such as merging of sources or rays (e.g. Pawlik 
Schay^|20Q8 ) and/or turning them off after a few Myrs. 


Turning them off seems acceptable, considering that the lu¬ 
minosity of a stellar population has dimmed by orders of 


magnitude 10 Myrs after its birth (see Fig. Al). However, 
Kannan et al. (2014b) have pointed out that the cumula¬ 


tive radiation from many such dim old sources may play a 
role in stellar feedback. Since we do not have an issue with 
the number of radiation sources in our implementation, stel¬ 
lar particles are never turned off after their birth, and the 
cumulative radiation from old populations is included. 

We implement three “separate” radiation feedback pro¬ 
cesses, describing different interactions between the radia¬ 
tion and gas: 

(i) Photons ionise and heat the gas they interact with, 
following the thermochemistry described in R13, typically 
heating the ionised gas to ~ 2 x 10^ K. All our runs with 
radiation include this process. We describe in Appe ndix [A| 
how we use the SED model from Bruzual & Chariot (2003) 


to derive photoionisation cross sections, which are typically 
a few times 10~^® cm^ for Hi, Hel, and Hell. 

(ii) Direct pressure, i.e. momentum transfer, from the ion¬ 
ising photons onto the gas. 

(iii) Indirect radiation pressure on the gas, via dust par¬ 
ticles, from the ionising photons, optical photons, and from 
reprocessed IR radiation, where the latter multi-scatters. 


R14 contains a detailed description of the implementation 
and tests of the latter two processes, including the diffusion, 
pressure, and work of multi-scattered IR radiation. We per¬ 
form the correct diffusion of IR radiation by a partition in 
every cell into free-streaming and trapped photons, where 
the trapped photons dominate in the case of large optical 
depth on the scale of the cell width. 

We will refer to the various radiation feedback processes 
under the collective acronym of RT (radiative transfer) feed¬ 
back. RT feedback may thus refer to the inclusion of any 


or all of the radiation feedback processes under considera¬ 
tion. We will successively add the three RT processes to our 
simulations, to probe their respective importance. Always 
included in RT feedback is photoionisation and photoionisa¬ 
tion heating, from the three UV photon groups. Onto that 
we add direct pressure from photoionisation (again only the 
three UV groups). Finally, we add radiation-dust interac¬ 
tions, from all five groups. 

Each photon group i has a dust-interaction opacity, ki^ 
listed in the rightmost column of Table The gas absorbs 
momentum from the photons (via dust) with the gas opacity 

sa^kiZjZQ, ( 6 ) 


i.e. in our model, the dust content simply scales with the 
metallicity of the gas. Higher energy photons (all but IR) 
absorbed by dust are reprocessed, i.e. re-emitted, into the 
IR group, while IR photons are (multi-) scattered by the 
dust. 

For the IR we assume an opacity of Km = 
10 Z/Zq cm^ , while for the higher energy photons we as¬ 
sume k\jy = 10^ ^/Z© cm^ i-e. a hundred times higher 
than that of the IR. These opacities are physically motivated 
from a combination of observations and dust-formation the¬ 
ory of the ISM and stellar nurseries ( Semenov et al.|2003 for 
IR, Li V Draine|2Q01 for higher energy radiation), but they 
are uncertain by a factor of few, due to model uncertainties 
and the temperature dependency, which we ignore. Similar 


values have been used in e.g. Hopkins et al. (2012c), Agertz 


et al. (2013), and Roskar et al. (2014). The IR opacity we 


use is at the high-end of what is usually considered in the 
literature, which is ^ir 5 — 10 cm^g“^. We have tested 
and confirmed that our results are insensitive to order-of- 
magnitude variations in the dust opacities (see ^4.3). 

There are two important exceptions from the default 
behaviour of the implementation described in R14. 

Firstly, our resolution of ~ 10 pc does not allow us to 
accurately capture the regime where dust is optically thick 
to photons, and radiation and dust are coupled via absorp¬ 
tion and blackbody emission. For this reason, and also for 
the sake of simplicity, we exclude the dust temperature evo¬ 
lution (§2.3.2 in R14), where the gas temperature is coupled 
directly to the IR radiation temperature via the Planck cross 
section. We decouple the dust temperature by simply setting 
the Planck cross section to zero (while keeping a nonzero 
Rosseland opacity). 

The second change is that we assume a fully directional 
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Figure 1. Illustration of radiation flux in the flve photon groups included in this work. The maps show density-weighted solid-angle 
integrated photon fluxes, cE^ along the LOS in the g9 galaxy with SN and full RT feedback (g9_SN_rhpd) at 250 Myr. The photon 
groups are shown by increasing photon energy, from left to right. The upper row shows the galaxy face-on and the lower row shows it 
edge on. The much larger contrast in the fluxes of the ionising photons (three right-most panels), owes to their much shorter mean free 
paths. Also, to a smaller degree, the optical photons have larger contrast than IR radiation, for the same reason. For the corresponding 
distribution of stars and gas in the same snapshots, this flgure can be compared to Fig. (bottom left panel). 


radiation in each cell for the free-streaming radiation pres¬ 
sure (R14, Eq. 28), by using a renormalised radiation flux 
magnitude of cE, rather than the actual radiation flux of 
|F| ^ cE. We do this to counter a resolution effect, as the 
reduced flux, = |F|/cE ^ 1, takes a few (~ 5) cell widths 
to evolve to unity with our advection scheme, even with free- 
streaming radiation. We demonstrate this numerical effect 
with a simple idealised experiment in Appendix [B] For the 
cell containing the emitting source, this resolution artefact 
is obvious, since the radiation is isotropic and hence has zero 
bulk flux (only E is incremented with stellar emission). The 
lack of bulk radiation flux very close to the emitting stellar 
particles diminishes the effect of radiation pressure, espe¬ 
cially since it turns out that HII regions are often poorly 
resolved in our simulations. Therefore, we apply this full 
reduced flux approximation = 1) for the radiation pres¬ 
sure, to compensate for resolution effects. It can then be 
argued that we overestimate radiation pressure, especially 
in regions where cancellation effects are relevant, but since 
it turns out that radiation pressure is very weak in our sim¬ 
ulations, we prefer to be in danger of overestimating rather 
than the opposite. We do not apply the full reduced flux 
approximation for the IR photon group, since pressure from 
the IR radiation, in the limit where the optical depth is not 
resolved, is accurately captured by the radiation trapping 
scheme (R14). 

Fig. □ illustrates the distribution of photons, for the 
five radiation groups, in one of our runs of the intermedi¬ 
ate mass galaxy disk (g 9). The figure shows mass-weighted 
averages along lines-of-sight (LOS) of photon fluxes, inte¬ 
grated over all solid angles, i.e. the mapped quantity is cE, 
where c is the reduced speed of light and E is the radia¬ 
tion energy density. From left to right, the maps show pho¬ 
ton groups with increasing energy, starting with IR on the 
far left, the optical, and finally the three ionising groups. 
The photon fluxes differ greatly between the photon groups, 
decreasing with increasing photon energy. We use differ¬ 
ent color scales, such that the logarithmic range is the 
same, but the upper limit roughly matches the maximum 
flux in each set of face-on/edge-on maps. For the high¬ 


est energy group (far right) the low luminosity is simply 
due to the low emissivity from the stellar populations (see 
Fig.|^ where we plot the emissivity of the stellar pop¬ 
ulations). For the two lower energy ionising groups (sec¬ 
ond and third from right) the stellar emissivity is similar 
to that of the optical group, yet the galaxy luminosity is 
clearly much lower than in the optical. This is due to the 
much more efhcient absorption of the ionising photons. For 
photoionisation of hydrogen and helium, the opacities are 
afrup ~ 6 X lO^cm^g”^, where a ~ 10“^® cm^ is the pho- 
toionisation cross section (see Table and Appendix [A| and 
nip is the proton mass, while for the optical group the opac¬ 
ity is hiopt = ^opt^/Z© = 10^ cm^ g“^. Hence the difference 
in opacities is more than three orders of magnitude. While 
the ionising photons are absorbed close to their emitting 
sources, the optical photons are much more free to propa¬ 
gate through the disk and escape from it. The direct stellar 
IR emission is relatively dim, about three orders of magni¬ 
tude lower than that of the optical group, yet the maps on 
the far left show that the radiation energy flux is highest 
in the IR group. This is because the IR photons are mostly 
reprocessed from the Optical and UV photons, which are 
captured by the dust and re-emitted into the IR. 

Fig. 0 illustrates the effect of the radiation on the hy¬ 
drogen ionisation fractions in the gas, which are tracked 
by the code. The left panel shows a run with SN feedback 
only, while the right panel shows the same galaxy also with 
full RT feedback, which results in an abundance of dense 
photoionisation-powered Hll regions. 

2.8 Overview 

Table ^ lists the properties of the simulated galaxies. We run 
each simulation for 500 Myr. Tablej^lists the 6 combinations 
of four feedback processes included in the simulations: No 
feedback at all (nofb), SN feedback only (sn), with added 
radiation feedback with radiation heating only (sn_rh), with 
added direct pressure from ionising photons (SN_RHP), with 
added radiation pressure on dust, and optical and (repro¬ 
cessed) IR radiation groups (sn_rhpd), and, finally, with all 
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Figure 2. Ionised hydrogen fractions in the g9 galaxy at 250 
Myr. The maps show mass-weighted ionised fractions along the 
LOS, for SN feedback only (left, g9_Sn) and added (full) radiation 
feedback (right, g9_SN_rhpd). The right map is the same snapshot 
as shown in the panels of Fig. [T] 


Table 3. Feedback processes included in the simulations 


Feedback 

SN 

Radiation 

Direct rad. 

Dust 

acronym 

feedback 

heating 

pressure 

pressure 

NOFB 





SN 

/ 




SN_RH 

/ 

/ 



SN_RHP 

/ 

/ 

/ 


SN_RHPD 

/ 

/ 

/ 

/ 

RHPD 


/ 

/ 

/ 


radiation feedback processes included, but without SN feed¬ 
back (rhpd). The name of each run is a combination of the 
acronyms from Tables and e.g. the name g9_SN_rhp 
represents the g9 galaxy (baryonic mass of 3.5 x 10^ M©), 
simulated with SN feedback and ionising stellar radiation 
with heating and direct pressure. 


3 RESULTS 

We now turn to the simulated galaxy disks and examine 
the impact of radiation feedback and the interplay of radia¬ 
tion and SN feedback. We start with the intermediate-mass 
galaxy, which has the highest resolution in terms of the num¬ 
ber of volume/particle elements, and then consider in turn, 
with somewhat less detail, the more massive GlO galaxy and 
the less massive g8 galaxy. 

3.1 G9: intermediate mass gas-rich galaxy 

We first focus on the intermediate mass galaxy, g9 (?^ one- 
tenth the baryonic mass of the MW), and we begin by con¬ 
sidering the qualitative effects of the different feedback pro¬ 
cesses on the morphology of the disk. Fig. shows maps 
of stellar density and total hydrogen column density, face 
on and edge on, at 250 Myr, which is half the run duration. 
Without feedback (top left panel), the galaxy contains many 


cold star-forming clumps interconnected by narrow gas fila¬ 
ments. SN feedback (top right panel) dramatically reduces 
star and clump formation, especially at large radii, smooths 
out the gas distribution and thickens the gas disk compared 
to the no feedback case. The inner ~ 3 kpc from the center 
of the galaxy remain quite clumpy, however. The addition of 
ionising radiation and photoionisation heating (middle left 
panel) adds to the effect of SN feedback by further smooth¬ 
ing the morphology of the galaxy, and further reducing the 
number of clumps. The addition of radiation pressure, direct 
(middle right) and on dust (bottom left), has little impact. 
With SN feedback excluded, radiation heating and pressure 
on its own (bottom right) is insufficient to prevent massive 
clump formation in the galaxy, and it is noticeably more 
clumpy and thinner than with SN feedback only. 

3.1.1 Star formation 

Star formation is the most direct probe of the efficiency 
of feedback processes. The more efficient the feedback, the 
more it will reduce and regulate star formation. Fig.j^shows, 
for the g9 galaxy, the cumulative stellar mass formed over 
time (upper panel) and star formation rates (solid lines in 
the lower panel). These results are in line with the qualita¬ 
tive effects we saw in the previous maps. Compared to the no 
feedback case, turning on SN feedback reduces the formation 
of stars by about 35% at 500 Myr. Turning instead to radia¬ 
tion feedback, with both the pressure terms included, gives a 
very similar reduction in the star formation. Combining SN 
and full radiation feedback (three thickest curves) consider¬ 
ably reduces the star formation again, by ~ 70% compared 
to the no feedback case, and by ~ 50% compared to the 
cases with SN or radiation feedback only. 

We can probe the importance of radiation pressure by 
comparing the curves where SN feedback is combined with 
successive introductions of radiation feedback processes, i.e. 
photoionisation heating, direct ionising radiation pressure, 
and radiation pressure on dust. The stellar mass formed is 
very similar, indicating that radiation heating is the domi¬ 
nant radiation feedback process. 

The < 10% variation in the stellar mass formed at the 
end of the runs for the various radiation feedback processes, 
is too slight to require serious interpretation. It is likely a 
random effect where small variations in the feedback model 
trigger massive clump formation at different times in the 
simulations. Individual clumps can dominate the star for¬ 
mation for tens of Myrs, while they migrate to the center 
of the disk. While these clump formations are likely ran¬ 
dom, we cannot rule out the possibility that these effects 
of successively added radiation feedback processes are sys¬ 
tematic. If the effect is real and non-stochastic, the way it 
works is somewhat counter-intuitive, as the addition of ra¬ 
diation pressure on dust and the sub-ionising photon groups 
(sn_rhpd) boosts star formation. This implies negative feed- 
baek, which can be explained by a scenario where radiation 
pressure sweeps the gas into concentrated star-forming shells 
or clumps. However, we do not see a negative feedback ef¬ 
fect from radiation pressure on dust in the other galaxies 
considered in this paper, and hence we conclude that it is a 
random effect, rather than systematic. 

Focusing on the star formation rates for the different 
feedback processes, in the bottom panel of Fig. we see 
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Figure 3. Maps of the g 9 galaxy (roughly ten times less massive in baryons than the MW) at 250 Myr, for the different feedback runs. 
Each panel shows face-on and edge-on views of the stellar density (left) and the hydrogen column density (right). From left to right, top 
to bottom, the panels show the runs without feedback (nofb), SN feedback only (sn), SN and radiation heating (sn_rh), SN-hradiation 
heating + direct pressure (sn_rhp), SN + radiation heating + direct + dust pressure (sn_rhpd), and radiation heating + direct + dust 
pressure (rhpd). The physical length scale and the color scales for the stellar and gas column densities are shown in the top left panel. 
The addition of radiation feedback smooths and thickens the disk, compared to SN feedback only. The respective additions of direct UV 
radiation pressure and then optical and IR pressure have little effect. 
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Figure 5. Outflows from the g 9 galaxy at 500 Myr. The maps 
show total hydrogen surface density for SN feedback only (left) 
and added (full) radiation feedback (right). The time, color-, and 
length-scales are marked in the left map. Dotted horizontal lines 
mark planes 0.2 -Rvir from the galaxy plane, where we measure 
the outflow/inflow rates plotted in Fig.[^ 


Figure 4. Star formation and outflows in the g 9 runs with dif¬ 
ferent feedback processes included, as indicated in the legend: no 
feedback (nofb), SN feedback only (sn), RT feedback only with 
all processes activated (sn_rh), SN+RT feedback with RT heat¬ 
ing only (SN_RHP), the same with added direct ionising radiation 
pressure (sn_rhpd), and then with added lower energy radiation 
and dust pressure (rhpd). Upper panel: stellar mass formed 
over time. Lower panel: star formation rates (solid lines) and 
outflows across planes at distances of 0.2 Rvir from the disk plane 
(dashed lines). 


that the combined SN+RT feedback flattens out the star 
formation history compared to the case of no feedback or 
individual SN or RT feedback. The star formation rates de¬ 
cline in the latter half of the runs with no feedback or in¬ 
dividual SN or RT feedback. This is due to the galaxy disk 
starting to be starved of gas. The initial disk gas mass is 
2 X 10^ M 0 , and it is clear from the upper panel that 
a considerable fraction of this mass has already been con¬ 
verted into stars at 500 Myr. This narrows the difference in 
star formation rates between the runs: while the rate is sup¬ 
pressed by a factor 2 at 500 Myr by the combination of 
SN and RT feedback, and not suppressed at all by only SN 
or RT feedback, the suppression factor is much higher before 
gas depletion sets in, peaking at a factor 5 at 150 Myr 
for the combined feedback case and a factor 2 — 3 for the 
“single” feedback case (excluding the first ~ 50 Myr, when 
the disc is relaxing). 


3.1.2 Outflows 


Galaxies produce outflows, and it has been suggested that 
radiation feedback, and radiation pressure in particular, may 


be important for generating these galactic winds (Murray 
et al.|[2QTT ). Fig. shows edge-on maps of the total hydro¬ 


gen column density for the g 9 galaxy at 500 Myr, with SN 
feedback only (left) and with added full RT feedback (right). 


The panels show that winds are generated in the g 9 galaxy. 
The winds are produced by SN feedback: maps (not shown) 
with no feedback or RT feedback only show no sign of winds. 
The figure reveals slightly different wind morphologies, with 
the SN_RHPD case showing a more collimated wind than the 
SN run, where the wind seems to form a conical shell, i.e. 
with a gas-free zone along the z-axis through the center of 
the disk. This difference is due to the star formation be¬ 
ing more concentrated towards the center of the disk in the 
SN+RT feedback case, while it is located in a few clumps at 
various radii from the center in the SN case. 

We consider the winds more quantitatively in the 
dashed curves in the lower panel of Fig. which show gas 
outflow rates across disk-parallel planes at |z| = 17.8 kpc, 
or 0.2Rvir, in each direction from the disk. The planes are 
indicated by dashed horizontal lines in Fig. These are 
gross outflow rates, i.e. we exclude from the calculation those 
cells intersecting the planes that have inflowing gas veloc¬ 
ity. Where outflows exist across those planes, which is in 
all the runs with SN feedback included, the outflow rates 
are similar, within roughly a factor of two of each other. RT 
feedback has very little effect on the outflow rates, regardless 
of whether or not radiation pressure is included. 

Outflows are often quantified in terms of the mass load¬ 
ing factor, which is the ratio between the outflow rate and 
the star formation rate. In the case of Fig.|^ the mass load¬ 
ing is quite low, i.e. the outflow rates are more than an order 
of magnitude less than the star formation rates. Although 
the outflow rates change little with the addition of radiation 
feedback, the mass loading is typically a few tens of percent 
higher, since the star formation is less efficient. 

In Fig. we focus on the end-time of 500 Myr and 
show gas flow rates and mean speeds across parallel planes 
as a function of distance |z| from the galaxy plane. Here 
we split the gas cells into those with outflowing and in¬ 
flowing z-velocities, shown in solid and dashed curves, re¬ 
spectively. RT feedback has very little effect on outfiow/in- 
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Figure 6. Outflow/inflow rates (upper plot) and speeds (lower 
plot) in the g 9 galaxy at 500 Myr. 


flow rates or speeds when added to SN feedback, except at 
1^1 < 3 X 10“^ i^vir ~ 300 pc, which is more or less inside 
the gas disk. At these small distances from the central plane 
of the disk, RT feedback slightly increases the outhow rates, 
but notably also the inhow rates, which follow the outhow 
rates closely. This shows that the RT feedback has the main 
effect of stirring up the disk gas without ejecting it from the 
galaxy. This matches with the qualitative differences in the 
edge-on maps in Fig. where the SN+RT feedback runs 
can be seen to have a slightly thicker and more diffuse disk 
than the SN only case. By itself, radiation feedback does 
not produce outhows (yellow curves in Fig. 1^, but it thick¬ 
ens the disk considerably compared to the no feedback case 
(green curves). 


3.1.3 The effect of the radiation 

We found in the previous subsection that radiation feedback 
helps regulate star formation in the g9 galaxy. Photoionisa¬ 
tion heating dominates the radiation feedback, while radia¬ 
tion pressure appears to have very little effect, if any. 

We now consider how the photons affect the proper¬ 
ties of the galactic gas. We compare in Fig. [^temperature- 
density phase diagrams of gas in the g9 galaxy, for the cases 
of no feedback (top left), SN feedback only (top right), RT 
feedback only (bottom left), and combined SN+RT feed¬ 
back (bottom right). For the RT feedback we have included 
ah radiation feedback processes, but we note that remov¬ 
ing radiation pressure, direct or on dust, has no discernible 
impact on the diagrams. 

The diagrams show stacked results from outputs every 


50 Myr for t = 100 — 500 My]0 starting after the initial re¬ 
laxation of star formation seen in Fig.^ We stack the results 
to show a crude time-average and reduce the stochastic inhu- 
ence of the formation and destruction of dense clouds, which 
can shift the maximum densities considerably. Apart from 
this shift in the maximum density tail, there is no qualita¬ 
tive change in the diagrams between the stacked snapshots. 
We will refer to results stacked by the same snapshots as 
“time-stacked” in the remainder of this paper. 

We over-plot the star formation thresholds in density 
and temperature (vertical and horizontal dotted lines), the 
median temperature per density bin (solid blue curve), and 
the mass-weighted mean density (solid blue vertical line). 
We bracket the mean densities by the maximum and min¬ 
imum means per stacked snapshot (dashed blue vertical 
lines), indicating the shift caused by the formation and de¬ 
struction of dense clouds. The diagonal dashed lines indi¬ 
cate the non-thermal Jeans pressure, Eq. ( [^, used to pre- 
vent resolution-induced fragmentation of gas ^Truelove et si. 


1997). The artihcial pressure term dominates the pressure 


of gas below this line, i.e. in the shaded region, making it 
the de-facto dominant “feedback” in this high-density low- 
temperature gas. Without feedback (top left), the Jeans 
pressure predominantly supports this gas, while adding SN 
feedback (top right), RT feedback (bottom left), or a com¬ 
bination of the two (bottom right), typically increases the 
temperature and decreases the density, and thus reduces 
the amount of gas supported by this artihcial pressure. One 
should keep in mind throughout that the effect of adding SN 
and RT feedback is somewhat weakened by the existence of 
this Jeans pressure, which must be in place in ah simula¬ 
tions as a last resort to keep gas from collapsing beyond the 
resolution limits. 

In the bottom left diagram, we see that radiation feed¬ 
back on its own increases the median temperature of dense 
gas compared to no feedback (top left), by heating a con¬ 
siderable amount of photo-ionised gas to ~ 10^ K. However, 
it has only a tiny impact on the mean density, compared to 
the no feedback case. Combined with SN feedback, radiation 
(bottom right) is much more efficient at decreasing the mean 
density, by almost half a dex compared to SN feedback only. 

Judging from the diagrams, the suppression in star for¬ 
mation due to radiation feedback appears to owe to either 
of two effects, or both: i) direct heating of the gas, which 
raises it above the temperature threshold of 3000 K for star 
formation, i.e. gas moves up, or ii) resistance to gas collapse, 
indirectly due to the heating, i.e. gas moves to the left. To 
investigate the direct effect, we have repeated runs g9_RHPD 
and g9_SN_rhpd, after removing the temperature threshold 
for star formation. The run with radiation feedback only, i.e. 
g9_RHPD, shows slight sensitivity to the temperature thresh¬ 
old, with 10% more stellar mass formed at 500 Myr with 
the threshold removed, while the run with SN+RT feedback 
(g9_SN_rhpd) actually produces 10% less stars if the tem¬ 
perature threshold is removed, which owes to an increase in 
the SN feedback efficiency. We conclude that the effect of ra¬ 
diation feedback is primarily due to adiabatic resistance to 


^ i.e. from outputs at t = 100,150, 200, 250, 300, 350, 400, 450, 500 
Myr 
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Figure 7. Temperature-density phase diagrams, time-stacked from snapshots every 50 Myr after relaxation, in the g 9 runs with different 
feedback processes, as indicated in the top right corner of each plot. The vertical solid lines show the mass-weighted mean density in 
each snapshot and the solid curves show the median temperature in each density bin. The dotted lines show the temperature threshold 
(horizontal) and density threshold (vertical) for star formation. The diagonal dashed lines indicate the non-thermal Jeans pressure, 
Eq. which is added to the other pressure terms (thermal pressure and trapped radiation pressure) in the hydrodynamics to prevent 
artificial fragmentation. Assuming a negligible contribution from trapped IR radiation, the pressure of gas below this line, as indicated 
by the shaded background color, is dominated by the artificial Jeans pressure, since the Jeans temperature is larger than T//r. 


gas collapse, rather than the precise temperature threshold 
for star formation. 

The bottom phase diagrams of Fig. [^reveal conspicuous 
features at the right end of the photo-ionised temperature 
plateau (?^ 2 x 10^ K), where the gas temperature decreases 
and density increases along narrow tracks. They are due to 
the HII regions being unresolved. The highest temperature 
tracks consist of single cells filled with radiation at a con¬ 
stant luminosity of a single young stellar particle, and can be 
accurately reproduced in single cell tests. The lower temper¬ 
ature tracks consist of cells adjacent to those source cells into 
which the constant luminosity propagates, again at roughly 
a constant rate. Under-resolved Hil regions are also visible 
in Fig.[^ indicated by a high contrast and “pixelated” peaks 
for the ionising photon groups in the three left-most panels. 
We will return to this resolution issue in §4.1[ 

Fig. shows the time-stacked mass-weighted density 
distribution of the gas in the g9 runs. Radiation and SN 
feedback clearly reduces the maximum gas density, but ra¬ 
diation pressure, when added, has very little effect. The plot 
supports the previous conclusion that the effect of the radi¬ 
ation heating is to prevent collapse of the gas by increased 



-4 -2 0 2 4 

log(nH [cm"^]) 


Figure 8. Time-stacked mass-weighted density distribution in 
the g 9 galaxy. Star forming gas is indicated by the shaded re¬ 
gion. SN and radiation feedback suppresses high gas densities. 
The suppression from radiation is dominated by radiation heat¬ 
ing, but IR and optical pressure on dust provides marginal extra 
support. 
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Figure 9. Comparison between the SN and SN+RT g 9 runs, of 
the probability distribution functions for the gas density at which 
stellar particles are created (upper panel) and produce SNe (lower 
panel). The shaded regions indicate star-forming densities. The 
solid curves in each panel show the cumulative probabilities. The 
upper panel indicates that RT feedback lowers the densities at 
which stellar particles are born, which should increase the SN 
feedback efficiency by allowing SN events to take place in a lower- 
density medium, as verified in the bottom panel. 


thermal pressure, which keeps the gas at lower densities. 
The effect is quite mild though, as the change in the density 
distribution is small when radiation is added to SN feedback. 


3 . 1.4 SN amplification 


Radiation can plausibly have the effect of amplifying SN 
feedback ( Pawlik Schay^|2009| Geen et al.||2015 ). It can 
diffuse the surrounding gas, which has the well-known ef¬ 
fect of decreasing the cooling rate, which scales with density 
squared. This in turn can make SNe more effective in stirring 
up the ISM, suppressing star formation and generating out¬ 
flows. This may happen as a combination of two effects: by 
preconditioning of the medium by the radiation before the 
SN events take place, but also in a preventive way, where 
the radiation feedback shifts the typical star formation den¬ 
sities to lower values, which directly causes SN events to 
take place at lower densities. 

Fig.j^shows the probability distribution of gas densities 
at which stellar particles are formed (upper panel) and at 
which they produce SNe 5 Myr later (lower panel), in the g 9 
runs with SN only and with full RT feedback added. From 
the upper panel we see that the RT feedback shifts star for¬ 
mation to lower densities, which now peak at the star forma¬ 
tion threshold, whereas they peak 1.5 dex above the thresh¬ 
old with SN feedback only. One also can read from the cumu¬ 
lative probability curves (solid lines) that with SN feedback 
only, about 45% of the stars are formed at nn ^ 10^ cm“^ 
(ten times the star formation threshold, n*), while 70% 


G9 



log(%,SN ! riH,b) 

Figure 10. Comparison between the g 9 runs with SN and 
SN+RT, of the probability distribution functions for the in¬ 
crease/decrease in the surrounding gas density between stellar 
particle birth and SN event. A value of zero at the x-axis in¬ 
dicates that the surrounding gas density stays unchanged from 
birth to SN, while negative/positive values correspond to a de¬ 
crease/increase in density. The solid curves show the cumulative 
probabilities. RT feedback has the effect of somewhat, but not 
dramatically, diffusing the gas around the stellar particles, before 
they produce SNe. 


of the stars form below the same density with RT feedback 
added. This then translates into a similar difference in the 
SN densities in the lower plot. With SN feedback only, about 
45% of the stars produce SNe in gas with densities below 
10 n*, while the addition of RT feedback increases this to 
70%. This similarity in the characteristic density difference 
between the two plots indicates that preconditioning of the 
gas by radiation does not play a major role. If it did, we 
should expect the typical SN densities to shift even further 
to lower densities. 

Even so, we take a closer look at the effect of radiation 
preconditioning in Fig. |10| Here we plot the probability dis¬ 
tribution functions, for SN and SN+RT feedback, for the 
relative difference between surrounding densities at stellar 
particle birth, nH,b, and SN event, nH,SN- The idea is that we 
remove the effect of the stars being born at lower densities 
with RT feedback. For the SN feedback only case, we find 
a strong peak in the probability around nu,SN = 1, 
which just means that typically a stellar particle’s birth and 
SN event happen at the same density. A slight majority, 
~ 60%, of the stars produce SNe at lower densities, and 
there is a tail in the distribution with a few percent of the 
SNe exploding at orders of magnitude lower densities. With 
radiation feedback added, the peak is still in the same place, 
but the distribution and the tail is shifted towards lower 
densities. The effect is not large though. 

In addition to giving information about the nature of 
(possible) SN amplification by radiation feedback. Figures 
and give us a hint about how radiation feedback sup¬ 
presses star formation. The radiation shifts star formation 
to substantially lower densities (Fig.|^, but does not as sub¬ 
stantially diffuse gas locally around young stellar particles 
(Fig, p^, suggesting that the effect of radiation feedback is 
more to prevent the formation of dense clumps, rather than 
destroying them after they form. 


© 0000 RAS, MNRAS 000 , 000-000 



































Galaxies that Shine 13 


Stars [M 0 pc‘^] 

Nh [cm-^] 

10 ' loVioV ‘ 10 -^- 10 “ 

» 

10 Kpc 250 Myr 

10 '’ 10 ” 10 ” 1 D”' 10 ” 10 ’“ 



G10_SN 



1 ♦ * 

■ -J ■' ■ i 

' ■ . 

♦ m ■ 

’ > ’ 

\ . 

A 



G10_SN_RHPD 



GIO 



time [Myr] 

Figure 12. Star formation and outflow rates in the GlO runs with 
different feedback processes included, as indicated in the legend. 
Upper panel: stellar mass formed over time. Lower panel: star 
formation rates (solid lines) and outffow rates across planes at 
distances of 0.2 i^vir from the disk plane (dashed lines). Feedback 
is much less effective here than in the less massive g 9 galaxy 
(Cf. Fig. 0. Radiation heating suppresses star formation more 
than SN feedback, but the effect is small. Radiation pressure is 
unimportant. Outffows are not affected by the radiation feedback. 


Figure 11. Maps of the GlO galaxy (baryonic mass of 3.5 x 
10^° Mq) at 250 Myr, for SN feedback only (upper panel) and 
full radiation feedback (lower panel). Each panel shows face-on 
and edge-on views of the stellar density (left) and total hydrogen 
column density (right). Radiation feedback has little noticeable 
effect in this galaxy, and in fact the same applies for SN feedback 
(the comparison to no feedback is not shown). 


3.2 GlO: Milky Way mass galaxy 

We now turn our attention towards our most massive galaxy, 
similar in mass to the Milky Way (MW) Galaxy. The galaxy 
is ten times more massive than the g9 galaxy we have anal¬ 
ysed so far, and of interest here is how the galaxy mass af¬ 
fects the radiation feedback. The mass is not the only thing 
different from the g9 galaxy, however. The metallicity of the 
gas is ten times higher and the gas fraction is considerably 
less: 30%, compared to 50% for the g9 galaxy. It makes sense 
to change also these properties, since the idea is to roughly 
follow the stages in the evolution of the present day MW. 
However, in §3.4| we will disentangle the effects of these dif¬ 
ferent galaxy properties on the radiation feedback. 

We first consider the qualitative effect of radiation feed¬ 
back on the galaxy morphology in Fig. EU where we compare 
face-on and edge-on maps at 250 Myr. We find no visible 
effect from radiation pressure, neither direct nor on dust, 
so we only compare here the case with SN feedback only 
(g10_Sn) and SN + full RT feedback (g10_SN_rhpd). The 


overall effect of adding RT feedback is less than in the g9 
galaxy (Fig. though the disk does become slightly less 
clumpy and more diffuse compared to SN feedback only. We 
do not show the no feedback case, but it looks very similar 
to the case with SN feedback, so also SN feedback is weak 
in this massive galaxy. 

We go on to compare the star formation, in Fig. |12| Here 
we again see that all the modelled feedback processes are 
much weaker than in the previous less massive galaxy. SN 
feedback initially slightly reduces the star formation com¬ 
pared to the no feedback case, but ends up with more stars 
formed (which is due to the recycling of gas in the SN case, 
resulting in an effectively larger gas reservoir). In such a mas¬ 
sive galaxy, SN feedback is though to become decreasingly 
important, and AGN feedback, which is not modelled, may 
start to dominate ( Bower et al.||2006 ). Also, it is likely that 
numerical overcooling becomes stronger, due to the increas¬ 
ing gravitational potential, gas densities, metallicity, and de¬ 
creasing physical resolution (although the larger stellar par¬ 
ticle mass should somewhat compensate by injecting more 
energy per SN event). 

In this galaxy, radiation feedback has a stronger effect 
on the star formation than SNe, though the effect is still 
weak, with an 7% reduction in the stellar mass formed (at 
500 Myr) with RT feedback only, and 10% if combined 
with SN feedback. The slightly increased feedback efhciency 
when combined with SNe hints at an amplihcation effect. 
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but the effect is small. Radiation heating dominates the ra¬ 
diation feedback, as adding radiation pressure and dust in¬ 
teractions has little effect on the star formation. 

Outflows rates across planes 0.2Rvir (38.4 kpc) from 
the disk are shown by dashed lines in the bottom panel of 
Fig. Outflows appear to be powered nearly exclusively 
by SN feedback, since the rates remain virtually unchanged 
after the addition of radiation feedback (of any sort). The 
mass loading factor of the outflow remains at < 0.1, similar 
to the g 9 galaxy. 


3.3 G8: gas-rich dwarf 


We now consider variations with RT feedback in the least 
massive galaxy, g8. Its properties only differ from those 
of the intermediate mass g 9 galaxy in terms of the halo 
and galaxy mass. The gas fraction and metallicity are un¬ 
changed, at 50% and 0.1 Z©, respectively. 

We begin with a qualitative comparison of morphologies 
with the inclusion of different feedback processes, shown in 
Fig. |13| We compare the cases of no feedback (top panel) full 
RT feedback (i.e. heating, direct and dust pressure, middle 
panel), and SN+RT feedback (bottom panel). RT feedback 
on its own is clearly much more efficient in this galaxy than 
in the previous, more massive ones. It completely suppresses 
the formation of massive clumps, smooths out density con¬ 
trasts, and considerably reduces the formation of stars. We 
do not show the case with SN feedback only, but note that 
in the weak gravitational potential of the g8 galaxy, it has a 
similar qualitative effect as RT feedback only, with the only 
clear difference being a somewhat thicker gas disk for SN 
only. Combining RT and SN feedback, however, has some ad¬ 
ditional impact on the galaxy morphology, with fewer stars 


and thicker, more diffuse gas disk (bottom panel of Fig. 13). 

We compare the star formation rates and outflows for 
the g8 galaxy in Fig. Here we see that the star forma¬ 
tion rates with RT feedback only are very similar to those in 
the SN only case. The combination of SN and RT feedback 
reduces the star formation by about 25% compared to in¬ 
cluding only one of those processes, which is much less than 
the relative reduction from the no feedback case when either 
process was added, which is 75%. In §3.1.3| we searched 
qualitatively for the existence of a feedback amplification in 
the g 9 galaxy, i.e. where the addition of one form of feed¬ 
back (RT) boosts the efficiency of another form (SNe) in 
quenching star formation, but found no clear evidence. Here 
we have an indication of the opposite effect. 

The inclusion of direct radiation pressure and dust in¬ 
teractions has no effect on the star formation rate. However, 
unlike the case of the more massive galaxies, it increases the 
mass outflow rates non-negligibly, restoring the outflow rate 
at late times back to that obtained with SN feedback only, 
as shown by dashed lines in the lower panel of Fig. |14| The 
effect comes predominantly from direct pressure from the 
ionising radiation, as can be seen by comparing the purple 
and dark red dashed curves. 


3.4 All galaxies: metallicity versus mass 

Comparison of the three galaxies in the previous subsections 
reveals a clear trend: the efficiency of RT feedback decreases 
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Figure 13. Face-on and edge-on maps of stellar density (left) and 
total hydrogen column density (right) for the g8 galaxy (baryonic 
mass of 3.5 X 10® Mq) at 250 Myr, with no feedback (top panel), 
full radiation feedback (middle panel), and added SN feedback 
(bottom panel). Radiation feedback alone efficiently prevents the 
formation of massive clumps. SN feedback alone (not shown) has 
a similar qualitative effect, though it results in a slightly thicker 
gas disk. Combining the RT and SN feedback smooths the gas 
distribution further. 
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Figure 14. Star formation and outflow rates in the g8 runs with 
different feedback processes included, as indicated in the legend. 
Upper panel: stellar mass formed over time. Lower panel: star 
formation rates (solid lines) and outffow rates across planes at 
distances of 0.2 i^vir from the disk plane (dashed lines). Radiation 
feedback is as effective as SN feedback at reducing star formation, 
but does it more smoothly, with SN feedback resulting in more 
bursty star formation. The relative effect of including a single 
feedback process (RT or SN) is much stronger than that of adding 
a second one. Radiation heating dominates the suppression in star 
formation, and reduces both the outffow rate and the mass loading 
factor of the outffow. Radiation pressure has a negligible effect on 
the star formation, but increases the outffow rate during the flnal 
250 Myr. 


with increasing galaxy mass. However, we varied not only 
the mass, but used a ten times higher metallicity in the GlO 
galaxy than in the less massive ones. We now want to inves¬ 
tigate how much the RT feedback efhciency is affected by 
galaxy mass, i.e. the gravitational potential, and how much 
by the gas metallicity, via its influence on the gas cooling 
time. For this purpose, we have run the three galaxies at both 
the metallicities we have considered, i.e. 1 Z© and 0.1 Z©. 

We quantify the efhciency of radiation feedback by cal¬ 
culating the relative reduction of stellar mass formed when 
a reference simulation is re-run with the addition of full RT 
feedback, i.e. 


c^(t) 


M.(t)x 


( 7 ) 


where M* (t) is the stellar mass formed in the simulation up 
to time t and X represents the feedback included in the ref¬ 
erence simulation. Values of < 1 correspond to feedback 
which suppresses star formation, while > 1 indicates neg¬ 
ative feedback, i.e. enhanced star formation. In Fig, we 
plot two such RT feedback efhciencies, upper 

panel, which shows the factor by which RT feedback sup¬ 
presses the stellar mass relative to the simulation without 



100 200 300 400 500 

time [Myr] 


Figure 15. RT feedback efficiency, i.e. cumulative suppression 
of star formation due to RT feedback (rhpd), plotted against 
time and compared for different galaxy masses and metallicities, 
as indicated in the legend. Different metallicities are denoted by 
solid and dashed curves, while galaxy mass is denoted by color 
and thickness, with increasing thickness indicating higher mass. 
An efficiency value of 1 corresponds to no effect on the star for¬ 
mation, while a value close to 0 indicates a strong reduction of 
the star formation. Upper panel: RT feedback efficiency when 
acting alone. Lower panel: RT feedback efficiency when com¬ 
bined with SN feedback. Both increased mass and metallicity re¬ 
duce the efficiency of radiation feedback (and SN feedback, not 
shown), except if radiation is combined with SN feedback, where 
the efficiency peaks for the intermediate galaxy mass. 


feedback, and in the lower panel, which shows the sup¬ 
pression when RT feedback is added to SN feedback. 

The upper panel shows the effect of radiation feedback 
in isolation, and gives a “cleaner” indication of the feedback 
efficiency than the lower panel, where the curves are quite 
sensitive to SN feedback efficiency, which is also (and in¬ 
dependently) sensitive to the galaxy mass, metallicity, and 
stellar particle mas^ However, the lower panel is quite im¬ 
portant, since the addition of radiation to SN feedback is 
more physically relevant than considering radiation feedback 
in isolation. We see from both panels that both increasing 
galaxy mass and metallicity weaken the effect of RT feed- 


^ Another factor, which we have not considered so far, is the effect 
that the stellar particle mass has on RT feedback. We have inves¬ 
tigated this for one of our galaxies, as discussed in Appendix [C] 
The indication there is that while stellar particle masses have a 
large effect on the SN feedback efficiency, they have much less 
impact on the RT feedback, which is likely because the energy 
injection is smooth rather than instantaneous. 
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Figure 16. Time-stacked mass-weighted probability distribution 
of IR optical depth, tjk along LOSs through the face of the GlO 
galaxy. The different curves represent inclusion of different feed¬ 
back processes, as indicated in the legend. A fraction of the gas 
in this galaxy is in the optically thick regime, tir > 1, where 
multi-scattering starts to play a role. However, it has little im¬ 
pact on the galaxy, as can be deduced from the similarity of the 
optical depth distributions with IR multi-scattering pressure ex¬ 
cluded and included (purple dot-long-dashed and solid dark red, 
respectively). 


back (the same applies for SN feedback, though this is not 
shown in these plots). 

The emerging qualitative picture is as follows: star for¬ 
mation in low mass galaxies is easily regulated by SN feed¬ 
back, due to a combination of long cooling times (low metal- 
licity and density), weak gravitational potential, and rela¬ 
tively massive stellar particles. Although RT feedback is by 
itself roughly as effective at regulating star formation (see 
Fig. 14), adding it to SN feedback has relatively little effect 
on the SF regulation, re duci ng the star formation rate by a 
few tens of percent (Fig. 15 lower panel)|^ With increasing 
galaxy mass, both SN and RT feedback become less effec¬ 
tive (Fig.j^and upper panel of Fig. |15| , but combining them 
may have a larger effect (Fig.[^ lower panel), though this 
is quite sensitive to the metallicity. At even higher mass, 
the gravity and cooling has become strong enough that not 
even the combination of feedback processes can significantly 
halt the star formation (Figs, and 15), especially at the 
higher, and more realistic, metallicities. 


3.5 IR multi-scattering 

The pressure due to multi-scattering IR photons has been 


cited as an important radiation feedback process ( 

Thomp- 

son et aL||2005| [Murray et aL[[2010[ [Hopkins et aL[ 

2012b|c[ 

Agertz et al. 

2013), yet we do not appear do get much of 


an effect at all from radiation pressure, including that of 


the IR radiation. In Fig. 16 we show the mass- (or column 
density-) weighted distributions of LOS IR optical depths, 
tir, through the face of the GlO galaxy, which has the largest 
optical depths of our disks. The LOSs are taken from time- 
stacked snapshots (every 50 Myr starting at 100 Myr, as 


[Hopkins et al.| ( |2012b| find qualitatively similar non-linear ef¬ 
fects when combining feedback processes. 


usual), and each has the width of the finest AMR reso¬ 
lution, or 36 pc. The plot quantifies how the mass is dis¬ 
tributed in face-on optical depths, which can safely be as¬ 
sumed to be consistently lower than edge-on optical depths, 
and thus more relevant for estimating the number of scat¬ 
terings (?^ Tir) photons typically experience before escaping 
the disk. 

A non-negligible fraction of the gas mass has larger than 
unity IR optical depths, so radiation trapping and multi¬ 
scattering does take place in the GlO galaxy, with maximum 
values of tir ~ 10. However, as we can clearly see by com¬ 
paring the curves in Fig. these opacities are not large 
enough for the IR radiation to diffuse the gas (and hence 
decrease the optical depths). Due to the lack of resolution, 
the gas does not reach the high densities, and hence optical 
depths, where multi-scattering plays a significant role. We 
can contrast these results to those of Hopkins et al. (2011), 
where typical optical depths around young stellar particles 
are found to be much higher, ~ 10 — 100, in HD simulations 
with ~pc resolution. 


Murray et al. (2011 ) argued that the collective radiation 
pressure from star formation can generate cold (~ 10^ K) 
outflows. Although our resolution is insufficient to resolve 
each individual optically thick cloud, we should in principle 
see this collective large-scale effect in our simulations, but we 
do not. The Murray et al. (2011) argument applies to mas¬ 
sive starbursting galaxies, and a critical star formation rate 
surface density of ~ 0.1 Mq yr“^ kpc“^. Our galaxies 
do reach E^c ~ 10 but this is confined to clumps and 

centers, with most of the disk below the critical value. How¬ 
ever, even if the star formation is mostly below the critical 
value, we would expect to see some effect of the radiation 
on outflows, and it is thus interesting that we see no clear 
effect at all. 

Resolution may still be the defining issue though: [Mur-j 


ray et al. (2011) envision neutral clouds where the radiative 


force acts on the surface facing the disc. We do not resolve 
these dense clouds, and radiation momentum is deposited 
more smoothly throughout whatever neutral gas exists in 
the halo. In §4.3| we explore qualitatively what we can ex¬ 
pect with better resolution, by artificially increasing the IR 
opacity (but find that outflows are still not generated). 


4 DISCUSSION 

Summarising the results, we find that radiation feedback 
has a modest effect on the star formation rates of our sim¬ 
ulated galaxies, while outflows are more or less unaffected. 
The suppression of star formation is due to the suppression 
of the formation of dense clumps. Radiation feedback be¬ 
comes less efficient with higher galaxy mass or metallicity, 
while the combination of radiation and SN feedback appears 
most effective at intermediate masses (and low densities). 
Photoionisation heating dominates the effect from radiation 
feedback, while radiation pressure, whether direct or from 
reprocessed, multi-scattering, IR radiation, has a negligible 
effect. 

We will now discuss several aspects of our findings, 
starting with the apparent inability to resolve Hil regions, as 
implied by Fig. We will then validate our results analyt¬ 
ically, comparing the relative impact of the different radia- 
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Figure 17. The plot shows the ratio of the Stromgren radius, 
rg, of ionised regions versus the maximum cell resolution, Ax = 
18 pc in g8/g9 (dotted blue) and Ax = 36 pc in GlO (solid 
red). The Hll regions are not resolved above the star formation 
density threshold (tt-h > 10 cm“^) in our simulations. The dashed 
green curve shows the ratio, at location of birth, if the total gas 
mass of a cell is always converted directly into a stellar particle 
(ignoring the mass depletion of the cell). It demonstrates that 
with the current star formation method, it is impossible to resolve 
Hll regions above the density threshold for star formation, if this 
treshold is > 50 cm“^. 


tion feedback processes, as expected in the numerical frame¬ 
work. Next, we will consider the expected effect of efficient 
IR feedback on star formation and outflows by artificially 
increasing the IR opacity. Finally, we will qualitatively com¬ 
pare our results to previous work on radiation feedback on 
galaxy scales, where the radiation effect is usually (but not 
always) modelled with subgrid recipes in pure HD simula¬ 
tions. 


4.1 On unresolvable Hii regions 


Our simulations show indications that Hll regions are not re - 
solved at gas densities nn ^ 10cm“^ (see Fig. [^and j3.1.3| l, 
which potentially affects our results at these high densities. 
We consider here in detail at what limit Hll region resolution 
becomes an issue. 

We can investigate this using the analytic expression for 
the Stromgren radius of a photo-ionised region in a uniform 
medium. The specific ionising luminosity of stellar sources 
is £uv ~ 5 X 10^® ionising photons per second per Solar 
mass, according to the Bruzual Chariot] (2003) model and 
assuming the |Chabrier| (2003) IMF (see Fig. Al). A stellar 
particle of mass m* then has luminosity Luv = 

(photons per second). The Stromgren radius around a stellar 
source is ( |Stromgren||1939| 

/ . \ 1/3 

_ / 3Luv 


rs 


= 21 pc 


(8) 


Cjjv 


5 X 1046 s-1 600 M© 


1/3 


( ) 

V10 cm“6 J 


-2/3 


10 cm“6, 

is the case B recombination 


where = 2.6 x 10 cm^ 
rate of hydrogen around 10^ K (Ferland et al. 1992), and 


where we have substituted the stellar particle mass used in 
the g8 and g 9 simulations, along with the star formation 
density threshold used in all our simulations. Fig, shows 
the ratio of the Stromgren radius and the cell width for all 
three simulated galaxies. For all simulations, the Hll regions 
around young stars are only resolved in gas below the star 
formation density, n*, and the stars must travel to densities 
of nu < 1 cm“^ within their “luminous” lifetime of ~ 5 
Myr to have their Hii regions resolved by more than ten cell 
widths. 

Judging from Eq. (|^, a simple solution to forcing Hil 
regions to be resolved in AMR simulations would be to in¬ 
crease the mass (and hence luminosity) of the stellar parti¬ 
cles. However, unless changes are made to the star formation 
recipe, we hit a concrete upper limit in which is the 

total gas mass of the hosting cell. 


m 


max 

* 


nurrip 

X 


{Axf 


1875 M© 


nn 

10 cm“^ 



where X = 0.76 is the hydrogen mass fraction. Assuming, 
for arguments sake, that we always convert the total gas 
mass of cells into stars, and that their neighbourhood has a 
roughly homogeneous gas density, the ratio of the Stromgren 
radius and cell width around newly formed stellar particles 
would be (from Eqs.[^and|^: 


rs 

Ax 


= 1.0 ( \ 

V50cm-3y' 


-1/3 


( 10 ) 


which is shown by the green dashed curve in Eig. Even 
at full gas conversion into stars, Hil regions are not well 
resolved for nn > 10 cm“^, regardless of the resolution. 

So far we assumed that the Stromgren sphere is powered 
by a single stellar particle. The situation changes when it 
becomes likely to have multiple young (^5 Myr) particles 
in a single cell, increasing the source luminosity and the 
size of the Hii region. Erom Eq. ([^, we can derive the star 
formation rate of a cell: 


SFRax = 1.75 X 10"® Mo/yr 


(V 

/ cff W riH y/2 

\18 pc/ 

h. 02 / vio cm- 3 / 


from which we can then derive the hydrogen number density 
at which more than one stellar particle, on average, is formed 
over 5 Myr, the time during which the stellar particles are 
luminous. This density, at which we can start expecting mul¬ 
tiple young stellar particles per cell, is 


mult -1 r‘^ 7 ■ —3 

riH = 167 cm 


600M© 



( 12 ) 


which coincides roughly with the density at which the track 
of temperature versus density widens to the right, in the bot¬ 
tom phase diagrams of Eig.[^ However, the phase diagrams 
demonstrate that, at the current resolution, the presence of 
multiple stellar particles in a single cell is insufficient to re¬ 
solve Stromgren spheres at high densities, as the cells with 
multiple stellar particles clearly do not reach T 10^ K. 

We now see that Hii regions cannot be fully resolved 
above these moderately large gas densities, unless changes 
are made to the star formation recipe, where e.g. more mas¬ 
sive stellar particles are formed from gas in a group of neigh¬ 
bouring cells or they are allowed to accrete gas in their life- 


© 0000 RAS, MNRAS 000 , 000-000 


































18 Rosdahl, Schaye, Teyssier & Agertz 


time. Stellar particle accretion is usually applied in numer¬ 
ical simulations of protostar formation and the evolution of 
individual molecular clouds, i.e. simulations at sub-galactic 
scales, and perhaps we are approaching a level of detail 
which requires some merging of methods for these different 
scales of galaxy evolution. Alternatively, one could apply 
stochastic radiation feedback, by allowing on average one in 
X particles to emit radiation at X times the default lumi¬ 
nosity. Such an approach has been used by |Dalla Vecchia fc] 


Schaye (2012) and Roskar et al. (2014) for SN feedback to 


overcome a related resolution problem of overcooling and an 
over-smooth distribution of stars. Stochastic radiation and 
SN feedback would thus appear to mesh quite naturally to¬ 
gether to overcome resolution problems. This is beyond the 
scope of the present paper though, and we can merely note 
the limitations in our feedback at high densities, which are 
in any case close to the resolution limit, where the pressure 
becomes dominated by the Jeans resolving pressure floor. It 
is presently unclear what the exact effect of under-resolved 
Hll regions is in our simulations, but likely it leads to an 
underestimate of the effect of photoionisation heating, since 
the gas in under-resolved Hll regions is heated to an unre¬ 
alistically low temperature (a fraction of the photo-ionised 
temperature which corresponds roughly to the ratio of the 
size of the real Hll region and the cell size). 


4.2 Analytic comparison of the RT feedback 
processes 

Among the main results of our simulations is that photoion¬ 
isation heating has a modest effect on regulating star for¬ 
mation, while radiation pressure contributes negligibly. We 
now seek to understand these results analytically, in order 
to see if they make physical sense and to ensure that they 
are not the product of implementation bugs. 

We can compare, within our numerical framework, the 
efficiencies of the different radiation feedback processes, 
i.e. photoionisation heating, direct pressure from ionising 
photons, direct pressure from optical photons, and multi¬ 
scattering pressure from IR photons. To simplify and quan¬ 
tify this comparison, we consider feedback in a single cell 
containing a radiation source and ignore radiation entering 
the cell from the outside. While we will write the follow¬ 
ing equations in terms of the simulation cell size. Ax, and 
the mass of stellar particles, m*, most of the equations also 
hold approximately for gas at a distance ~ Ax from a star 
(cluster) of mass m* with the assumed (and theoretically 
motivated) specific luminosity, provided the density, tem¬ 
perature, and metallicity are nearly uniform within Ax. 

We compare the radiation feedback efficiencies in terms 
of approximate “effective” temperatures. For photoionisa¬ 
tion heating, this is equal to the temperature of gas pho- 
toionised by stars, while for direct and IR radiation pres¬ 
sure (with the motivation of comparing those processes in a 
simple way), it is dehned by equating the radiation pressure 
and the thermal pressure. 


extends outside the cell, then the cell is simply heated to 
Thii, otherwise it is heated to a fraction of that tempera¬ 
ture which reflects the ratio of the volume of the Stromgren 
sphere to that of the cell, i.e. the host cell is heated to 


' 2 X 10^ K X min (/voi, 1), 


(13) 


where 


/vol = 


(Ax)^ 


= 6.7 


jCuv 


m. 


5 X lO^e s-iMq^ 600 Mq 

-3 


(14) 


( 1 

rV ^ 

V10 cm“^ / 

vispcy 


and we substituted our (g8 and g 9) simulation parameters 
for Ax, £uv, and m* (and the star formation threshold 
for uh). The specific stellar population luminosity, for the 
UV and for the other photon groups, can be read (approxi¬ 


mately) from Fig. Al 


4-2.2 Direct pressure from photoionisation 


To quantify the effect of radiation pressure and compare 
it to photoionisation heating, we measure it in terms of an 
effective temperature, corresponding to the pressure applied 
via momentum absorption from the radiation, and dehned 
as 

Teb = (15) 

pKB 

The radiation pressure is roughly the momentum absorption 
rate in the cell, p (momentum per unit time), divided by the 
cell area. 


Rrad = 


P 

6 (Ax)^ 


(16) 


The momentum absorption rate can be estimated from the 
luminosity of the stars contained in the cell and the opacity 
of the cell gas. The effective temperature is approximate, 
because we neglect the dependence on the mean gas particle 
mass /i, and we assume the radiation pressure to be isotropic. 

We assume for simplicity that the cell gas is in pho¬ 
toionisation equilibrium with the emitted radiatior0 and we 
ignore the radial dependence of the neutral fraction inside 
resolved Hll regions. We can then use the size of the pre¬ 
dicted Hll region, given by Eq. (|^, to estimate the fraction 
of the ionising luminosity contributing to the direct radia¬ 
tion pressure in the emitting cell, giving 


rrUV 


Ljjv 


c 6{Axf pkB 


1.2 X 10'^ K 


( rV 

V10 cm“^ / \ 


xmin(/^J, 1) 
Djjv 


(17) 




2 X 10^® ergs“^ ^0 

-2 

X min (/;;}, 1), 


Ax 
18 pc 


where we now measure the luminosity (and specific), Luv 
(>Cuv), in terms of energy rather than photon count (the 


4-2.1 Photoionisation heating 

Photoionisation tends to heat the ionised gas to Thu ~ 
2 X 10^ K, as we have seen in Fig.[^ If the Stromgren radius 


^ We ignore the instantaneous pressure from the radiation when 
the cell is in the process of being ionised, leading us to under¬ 
estimate the direct pressure at low densities, where growing Hll 
regions are resolved. 
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value is again typical for the SED model utilised). The right¬ 
most term accounts for whether the Hll region is resolved 
or not: the fraction of the ionising luminosity pressurising 
the emitting cell is the volume of the Hii region over that of 
the cell, but this fraction is roofed at unity, meaning all the 
emitted photons are absorbed in the emitting cell as the Hll 
region becomes unresolved. 


4-2.3 Pressure on dust from optieal photons 


The effective temperature corresponding to the pressure on 
gas via dust from optical photons is 


rjiOpt 


Lppt 1 rup /-I _ -ropt^ 
C 6 (Ax)^ pkB 


(18) 


'l.SxlO^K -- T—^ 

3 X 1036 ergs-i 600 ^ 


X 


( rb 

^ Ax \ 

V10 cm“3 J ' 

\18 pc) 


where ropt is the optical depth of the host cell: 


Topt ~ /^optpAx (19) 

^12 ^Opt ^ riH Ax 

10^ cm^g-i Zq 10 cm“^ 18 pc’ 

We ignore pressure on dust from UV photons, because 
the pressure from the UV photons is already counted, in 
for photoionisation, for which the opacity is orders of 
magnitude higher than for dust absorption. 


4 . 2.4 Multi-seattering pressure on dust from IR photons 


The effective temperature for multi-scattering reprocessed 
IR photons is 

Lir 1 TTip _ Lir Kmmp 


rpIR- 


c 6 (Ax) pkB 


tir = 


22 K 


-^Opt 


c GAx/cb 


( 20 ) 


3 X 1036 ergs-1 OOOM© 
^iR Z f Ax 


10 cm^ g“i Z 


o 


18 pc 


We only consider the optical stellar luminosity, since the IR 


luminosity is negligible in comparison (see Fig. A1). The IR 


multi-scattering feedback depends on the optical photons 
being absorbed and re-emitted into the IR. It is a safe as¬ 
sumption though, that this is true under any circumstances 
where multi-scattering is important, since Km ^ siOpt- 

The expression for the IR effective temperature assumes 
trapping of photons originating within the cell and ignores 
additional trapping of photons originating from the neigh¬ 
bouring environment. The previous expression should there¬ 
fore be taken as a lower limit. This is less of a concern for 
the other radiation feedback processes, since they are unre¬ 
solved at high (star-forming) densities, and we thus expect 
much less inter-cell flux of photons. 


4 . 2.5 Relative impaet of the radiation feedbaek proeesses 

To compare the radiation feedback processes (Eqs. |13|20| |, 
we replace the stellar mass, m*, by the fraction ^ 1 of the 
gas mass in a cell at a given density and volume. The value 
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/* = 1 gives an approximate upper limit on the effective 
temperature estimates for each of the processes, with the 
possible exception of the IR radiation, where things are more 
uncertain, as argued in the previous sub-section. We thus 
assume 




nnrup 

X 


{Axf = 1875 Mo f. 


hh 


10 cm“3 


/ Ax 
\ I8pc 


3 


( 21 ) 


In practice, the mass of each stellar particle is some fraction 
(1/3 in the g 8 and g 9 simulations, 2/3 in GlO) of the cell 
gas mass at the density threshold for star formation. The 
upper limit we use in Eq. (21), reflects the fact that at large 
densities, multiple stellar particles can be expected to form 
in the same cell over a short timescale, but no more than 
the total gas mass in a cell can be converted into stars. 


Substituting Eq. (21) into (I3)-(20) gives: 


= 2 X 10^ 


TTA = 3.7 X 10^ K 


K X min(/voi, I), 
Cjjv 


Ax 


(/> 


vol 5 


2 X 1036 ergs“^ 

1 ) /*, 


= 5.6 X 10^ K 


-^Opt 


Ax 


(1-, 


= 69 K 


3 X 1036 ergs“3 Mq^ 18 pc 
l^Opt Kir Z 


3 X 1036 ergs“^ Mq 


10 cuP g“i Z(7 


© 


riB 


Ax 


I0cm“3 yis pc 
where the volume fraction. 




( 22 ) 

(23) 

(24) 

(25) 


/vol = 


Tin 


2.1 X 102 cm- 


jCuv 


5 X 10-^6 s-im;; 


(26) 


is independent of resolution. 

These effective temperatures are plotted in Fig. for 
/, = 1 and Z = Zq, with the thick curves representing 
the 18 pc resolution used for our lower-mass galaxies, and 
the thin curves corresponding to Ax = 1 pc. We also plot 
(in dotted grey) the artiflcial non-thermal pressure, Eq. §, 
that is imposed to resolve the Jeans scales, and we shift it 
with resolution according to Eq. (|^. 

Photoionisation heating (solid red) dominates at low 
densities, nn ^ I0^cm“^, heating the gas to ~ 2 x 10^ K, 
while for higher densities, radiation feedback is surpassed 
by the Jeans pressure. In the absence of the Jeans pres¬ 
sure, IR multi-scattering would dominate at high densities. 
The direct UV and Optical radiation effective temperatures 
plateau at high densities, as the total particle luminosity 
becomes absorbed in the local cell, and increasing the reso¬ 
lution only makes radiation pressure weaker, since the lower 
stellar mass (and hence luminosity) has a stronger negative 
effect ((X (Ax)^) compared to the positive effect of the de¬ 
creased cell area (oc (Ax)“^). The IR radiation pressure is 
the only term which keeps rising for increasing densities, 
which is due to the multi-scattering, and it dominates over 
other radiation feedback processes at extreme densities. The 
IR effective temperature is, however, lower than the artifi¬ 
cial pressure floor, Tj, by about an order of magnitude (at 
Z = Z©), at any resolution. 































20 


Rosdahl, Schaye, Teyssier & Agertz 



Figure 18. Effective temperatures acquired iu a siugle emittiug 
cell via the differeut radiatiou feedback processes (Eqs. |22|25| >, 
assumiug Z = Zq aud that a fractiou /* = 1 of the cell gas is cou- 
verted iustautaueously iuto stars. For the thick hues, we assume 
Ax = 18 pc, while for the thiu hues we assume 1 pc resolutiou. 
The processes cousidered are: photoiouisatiou heatiug (solid red), 
direct pressure from photoiouisatiou (dashed blue), direct pres¬ 
sure via dust from optical photous (dot-dashed greeu), aud multi- 
scatteriug pressure from IR photous (dot-dot-dot-dashed yellow). 
Also plotted, iu dotted grey, is the effective temperature of the 
deusity-depeudeut pressure floor (Eq.|^. Photoiouisatiou heatiug 
domiuates the radiatiou feedback at low deusities iu this siugle 
cell limit, while the pressure floor takes over at high deusities. 
The first-order effect of iucreasiug the resolutiou is to decrease 
the effect of radiatiou pressure. However, at extreme deusities, IR 
trappiug by multiple stellar particles ou scales larger thau the cell 
width is likely to give a boost over that iudicated iu the plot. 


Fig. [^qualitatively justifies the results of our simula¬ 
tions. Comparing the plot to the bottom right phase diagram 
of Fig. [^ the plateau of 2 x 10^ K gas is clear in both 
figures, and the drop-off in temperature, which is a manifes¬ 
tation of unresolved Hll regions, occurs at a similar density, 
though slightly lower in the phase diagram, which is because 
less than the full cell mass is converted into stars in the simu¬ 
lations (/hc = 1/3). Radiation (heating) feedback is effective 
in preventing gas at low densities from clumping, but futile 
in dispersing clouds once the densities become high. Radi¬ 
ation pressure vanishes at low densities, and is negligible 
compared to the artihcial Jeans pressure at high densities. 

Judging from Fig. it appears that the modelled ra¬ 
diation pressure is doomed to always remain weaker than 
the artihcial pressure hoor we are forced to apply, especially 
considering that we are assuming an extreme upper limit 
where the full gas mass in a cell is converted instantaneously 
into stars. However, we stress again that we only consider in 
this analysis the effect in a single cell, and ignore the effect 
appearing from external stellar populations in neighbouring 
cells (and we also ignore the fact that particles can move 
to higher or lower densities during their lifetime). We can 
conclude that direct radiation pressure is weak at any res¬ 
olution, but with many stellar particles forming in highly 
resolved optically thick regions, we may see a considerable 


boost in the pressure from trapped multi-scattering IR ra¬ 
diation with higher resolution (and only at high densities). 
It remains a task for future work to establish what kind of 
resolution is required to see such a boost, but in the next 
subsection we will investigate what effects we can expect 
from it on large scales. 

The above analysis does not apply to the time- 
integrated effect of collective, direct long-range radiation 
pressure from many stellar populations on galactic scales, 
e.g. in stirring diffuse gas or pushing cold clouds out of the 
galaxy ( Murray et al.|20lT ). However, this effect is unimpor¬ 
tant in our simulated galaxies, since it exists (and is in fact 
exaggerated due to our full reduced hux approximation), yet 
radiation feedback is dominated by heating, with radiation 
pressure at best having a marginal impact. 


4.3 What does it take for IR radiation pressure to 
dominate, and what happens then? 

Up to this point, we have found that IR radiation has only a 
marginal effect on our galaxies, and we have shown analyti¬ 
cally that these results are to be expected with our current 
model and resolution. The considerable increase in resolu¬ 
tion, that appears to be required to investigate IR pres¬ 
sure feedback on small scales, and possible cascading effects 
on larger scales, is beyond our reach in the current paper. 
We can, however, instead artihciahy allow the IR radiation 
to dominate the galactic feedback by simply increasing the 
IR opacity. This can give us an estimate of how far we are 
from efhcient IR feedback, and, more importantly, how the 
galaxy reacts when IR feedback does become efhcient on 
small scales. Thus we get an idea about what to expect if 
we resolve the very large optical depths that are required 
for multi-scattering radiation pressure to play a role. For 
example, does the radiation generate large-scale winds, and 
does it create a much thicker gas or stellar disk? We thus 
ran variants of the gI0_RHPD simulation (where the IR op¬ 
tical depths are largest) with increased IR opacities. We 
compare here results for ^ir = (10^, 10^, 10^) cm^ g“^, i.e. 
ten, a hundred and a thousand times the default, physically 
motivated, opacity that we have used so far. We also set 
^uv = ^Opt = ^iR in the highest opacity run in order to 
increase the IR reproduction from higher-energy photons in 
line with the opacity increase. 

Fig. shows star formation rates and outhow rates 
across planes 0.2 Rvir from the disk plane. We hardly see 
any effect on the star formation rate, though the outhow 
rate is slightly reduced. Further increased opacity increas¬ 
ingly suppresses both star formation and outhow rates, with 
a very bursty star formation and almost totally quenched 
outhows in the most extreme case. The reduction in the 
outhow rate, which is at least as large as for the SFR, even 
with the appearance of bursty star formation, strengthens 
the impression of a non-violent radiation feedback, which 
stirs up the gas but does not systematically eject it. 

Fig. shows the ehect of the increased opacities on the 
gas density distribution. As expected, the IR radiation sup¬ 
presses high densities more efhciently with increased opacity. 
For the highest opacity, the density distribution cuts oh at 
the star-format ion density threshold (n* = 10 cm“^). This 
indicates that the IR pressure does indeed become very ef¬ 
hcient at preventing gas to form stars, but that is more or 


© 0000 RAS, MNRAS 000 , 000-000 









Galaxies that Shine 21 


G10 SN RHPD 



time [Myr] 


Figure 19. Star formation rates (solid lines) and outflow rates 
across planes at distances of 0.2 J?vir from the disk plane (dashed 
lines) in the g10_SN_rhpd galaxy with increased IR opacity. The 
thinnest (dark red) curve shows the results for the default opacity 
that we have used so far and the successively thicker curves show 
results where the IR opacity is increased, each time by a factor 
of ten. The star formation becomes bursty in the case with the 
highest opacity. Outflow rates decrease with increasing opacity, 
more or less in line with the reduced star formation, indicating 
that the radiation disrupts star-forming clouds gently, rather than 
violently. 


G10_SN_RHPD 



log(nH [cm ^]) 


Figure 20. Time-stacked mass-weighted gas density distribution 
in the g10_SN_rhpd galaxy, with increasing IR opacity. For the 
highest opacity used, when the star formation is reduced by more 
than an order of magnitude, the density distribution becomes cuts 
off at the star formation threshold, nu = I0cm“^. 


less the whole effect, i.e. the gas is kept diffuse, but within 
the galaxy. The IR pressure shifts SN explosions to lower 
densities, but this does not lead to increased outflow rates, 
as the reduced star formation more than compensates to 
reduce the outflows, leading to a decrease in mass loading 
with increasing opacity. 

Finally, we compare the galaxy morphologies for the dif¬ 
ferent opacity values in Fig. |21[ The galaxy simply becomes 
smoother with increasing opacity, both in the gas and stars. 





kjf,=10^ cm^/g 



^ "J 



kjR=10^ cm^/g 




Figure 21. Maps of the GlO galaxy (3.5 x 10^® Mq in baryons) at 
250 Myr, for SN and full RT feedback (same as bottom left panel 
in Fig. IF). but with increased IR opacity, as indicated in each 
panel: the opacities are increased from the default by a factor 10 
(top), 10^ (middle), and 10^ (bottom). The increased dominance 
of IR radiation pressure simply has the effect of smoothing out 
the peaks in gas density, and suppressing star formation. 


© 0000 RAS, MNRAS 000 , 000-000 












































22 


Rosdahl, Schaye, Teyssier & Agertz 

Comparison with other work 


4.4 


While many studies exist in the literature where radiation 
feedback (on galactic scales) is modelled with subgrid recipes 
in pure HD simulations, there are only a few in which a sub¬ 
set of the radiation feedback mechanisms that are modelled 


here are studied with RHD (jPetkova & Springel|2011 

jWise 

et al.||2012al |Kim et al.||20131 jHasegawa & Semelin 

to 

o 

CO 

Pawlik et al.|2013| 2015). 


In cosmological simulations of reionisation-era galaxy 
formation, Wise et ah (2012a) found direct radiation pres¬ 
sure to play a major role, suppressing star formation 
strongly and boosting outflow rates, but did not report on 
the isolated effect of photo-ionisation heating, which was al¬ 
ways included, and they did not include IR radiation effects 


(which are likely weak in such metal-poor galaxies). Petkova 
Spring'd (2011); Hasegawa & Semelin ( |2013| ; [Pawlik et al. 


(2013) also modelled reionisation-era galaxies, but only in¬ 
cluded radiation heating. They all found that radiation heat¬ 
ing gently suppressed star formation, while they did not re¬ 
port on boosted outflows. Pawlik et al. (2015) considered 
the addition of SN feedback and found that it dominated 
over the effect of radiation heating on star formation histo¬ 
ries. 


Kim et al. (2013) used the RHD implementation from 
Wise et al. ( 2012a| ) on an isolated ~ MW mass galaxy, and 


found slight (~ 20%) suppression of star formation, which 
they attributed to radiation heating, rather than radiation 
pressure. 


A large amount of work exists where radiation feedback 
has been included in pure HD simulations in the form of 
subgrid recipes, which are often empirically motivated. Al¬ 
though there are quantitative, and sometimes qualitative, 
differences, these studies broadly agree that IR radiation 
pressure strongly suppresses star formation and generates 


20T3l IHopkins et al.||2011| |2012c|b|a| 

Aumer et al.j 2013] 

Agertz et al.||2013| j Agertz & Kravtsov 

2015] Roskar et al.j 


2014). There are exceptions though: Ceverino et al. (2014) 


and Moody et al. (2014) found direct radiation pressure to 


mildly suppress star formation, while radiation heating and 
IR radiation pressure had a negligible effect. [TVujillo-Gomez] 
et al. (2015), on the other hand, found radiation heating 


dominated over direct radiation pressure in suppressing star 
formation, while outflows were not affected by the radiation 
(and IR effects were not considered). 


Our results do not show a wide and general agreement 
with previous studies of the effects of radiation feedback on 
galactic scales, which is not surprising, since there is no gen¬ 
eral agreement in the literature. The discrepancies probably 
largely come down to resolution. It appears that both RHD 
and HD simulations that show a substantial effect from IR 
radiation pressure have either sub-pc resolution or a subgrid 
model that boosts the optical depths. We lack sub-pc res¬ 
olution in the current paper, and we have so far made no 
attempt to compensate this with a subgrid model. Of these 
two options, we prefer in future work to increase the reso¬ 
lution, to probe from first principles how radiation feedback 
affects small scales, and how this effect may (or may not) 
cascade to larger scales. The strongest general disagreement 
we can find with other work concerns outflows. Where they 
are studied in the literature, radiation feedback appears to 
boost outflows most of the time, which is in contrast with 


our simulations. Our experiments with boosted IR radiation 
opacities hint that increased resolution will still leave us with 
a lack of radiation-generated outflows, but in the end, the 
best way to find out is to actually increase the resolution. 


5 CONCLUSIONS AND FUTURE WORK 


We ran and analysed adaptive mesh refinement simula¬ 
tions of isolated disk galaxies of baryonic masses 3.5 x 
(10®, 10^, 10^°) M© (the largest mass being comparable to 
that of the Milky Way), using a maximum resolution of 
18 pc. We studied the effects or stellar radiation feedback, 
which was modelled with radiation-hydrodynamics, acting 
on its own and also combined with (“thermal dump”) su¬ 
pernova feedback. We compared the effects of three separate 
radiation feedback processes: photoionisation heating, direct 
radiation pressure from UV and optical photons, and pres¬ 
sure from multi-scattered, reprocessed IR radiation. These 
are the first galaxy-scale simulations which model all these 
processes concurrently and with RHD. Our main findings 
are the following: 


• Stellar radiation feedback suppresses star formation in 
the simulated galaxies. It does so predominantly by pre¬ 
venting the formation of star-forming clumps, rather than 
by destroying those that form. The suppression of star for¬ 
mation with radiation feedback (ranging from a factor of 4 
for the low mass galaxy to only ~ 0.1 for the most massive 
one) is similar to that of “thermal dump” SN feedback. 

• Radiation feedback does not significantly amplify the 
efficiency of SN feedback, and in fact there is a hint of the 
opposite effect in the lowest-mass galaxy we consider, where 
the combination of radiation and SN feedback results in a 
weaker star formation suppression than one would naively 
expect from multiplying the individual suppression factors, 
although the combined effect does exceed that of the indi¬ 
vidual feedback processes. 

• Radiation feedback has a negligible effect on galaxy out¬ 
flows. If anything, the outflow rates are slightly suppressed, 
owing to the reduced star formation and subsequent decrease 
in SN activity. The outflow mass loading factor, i.e. the ra¬ 
tio between the outflow rate and the star formation rate, is 
typically of the order of 10~^, which is very low compared to 
non-RHD simulations that use subgrid recipes for radiation 
feedback (e.g. Hopkins et al.||2012b ). 

• As with (“thermal dump”) SN feedback, the effect of 
radiation feedback on star formation weakens with galaxy 
mass and metallicity. The combined effect of SN and radia¬ 
tion feedback is strongest in our intermediate-mass galaxy, 
which has a baryonic mass of 3.5 x lO^M©, i.e. about one- 
tenth of the mass of our Milky-Way. 

• The dominant form of radiation feedback is photoion¬ 
isation heating, while the effect of radiation pressure, both 
direct and on dust, is borderline negligible. We are able to 
explain the relative efficiencies of the different radiation feed¬ 
back processes using simple analytic estimates within the 
context of our numerical models. 

• The analytic estimates suggest that the effect of di¬ 
rect radiation pressure from ionising radiation on galaxies 
is likely negligible in reality. However, multi-scattering ra¬ 
diation pressure from IR radiation is not properly captured 
in our simulations. This is because our resolution (~ 10 pc) 
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does not allow the collapse to sufficiently large local densi¬ 
ties for the gas to become significantly optically thick to the 
IR radiation. 

• To estimate the qualitative effects of multi-scattering IR 
radiation that may be revealed by future, higher-resolution, 
simulations, we ran simulations in which the IR opacities 
were boosted by orders of magnitude over realistic, physi¬ 
cally motivated, values. This boost makes multi-scattering 
radiation effective at regulating star formation, but in a gen¬ 
tle way that merely smooths out the galaxy disk, without 
generating outflows. 

• Resolution is also an issue for the ionising radiation. 
With the current method for forming stars, where stellar 
population particles are instantaneously formed out of the 
gas in a single AMR cell, Hil regios are unresolved at den¬ 
sities riH ^ 10cm“^, regardless of the numerieal resolution. 
The consequence of the unresolved Hll regions is likely an 
underestimate of the regulating effect of radiation heating on 
star formation in dense gas, since the gas cells hosting young 
stars are heated to temperatures lower than the ionisation 
temperature. Possible ways to deal with this problem in the 
future include stochastic radiation feedback or a modified 
method for star formation. 

• Although we have not considered this in detail, we find 
in the resolution tests described in Appendix that radi¬ 
ation feedback is much less sensitive to the stellar particle 
mass than is (“thermal dump”) SN feedback. This makes 
sense, since the radiation is continuous, while the SNe are 
instantaneous, and for explosive feedback the radiative losses 
decrease for higher maximum temperatures. 


An important caveat for our study is that our simulated 
galaxies do not have the high surface densities that occur in 
the massive, starbursting galaxies that have been the focus 
of theoretical work which predicts efficient regulation of star 


formation and outflows by radiation pressure (e.g. Murray 
eTalllMT l. At high-redshift {z ~ 3), where gas accretion 


and star formation peak, radiation pressure may even play 
a role in ‘normal’ low-mass galaxies. In the future we will 
expand our simulations to include more massive galaxies, 
and gas-rich galaxies representative of high redshift, which 
may exhibit greater sensitivity to radiation pressure. 

We also note that the choice of SN feedback recipe likely 
affects the interplay of feedback processes and the net effect 
of radiation feedback. For simplicity, and in order to make 
sure we did not over-inject feedback energy in this first round 
of simulations, we used “thermal dump” SN feedback, which 
is known to be inefficient and suffer from resolution-induced 
overcooling. In future studies it will be interesting to see how 
the interplay of feedback processes is affected by the use of 
more efficient (and more realistic) SN feedback recipes. 

There are many additional interesting paths to follow, 
such as improvements of our radiation feedback model, the 
inclusion of other sources of radiation than stars, and an 
expansion to both larger and smaller physical scales. 

An interesting model improvement is to consider the ef¬ 
fect of the local radiation field on metal cooling, which has 
been suggested by Cantalupo (2010) to effectively quench 
cooling and subsequently star formation in galaxies. Another 
important model improvement is the inclusion of the forma¬ 
tion and radiative dissociation of H2, which is highly rel¬ 
evant for studying star formation in detail. AGN feedback 


may be fundamentally radiative in origin, and it is quite 
interesting to see in what ways, if any, RHD experiments 
would differ from subgrid recipes. It is relatively straightfor¬ 
ward to add AGN radiation to our simulations, as long as a 
recipe for black hole accretion is in place. Some additional 
radiation processes, such as Gompton scattering, are likely 
important, and it is quite likely that it will remain difficult 
to resolve optically thick regions properly. 

We intend to study radiation feedback on scales both 
larger and smaller than the current study. The larger scales 
involve cosmological zoom RHD simulations, where the ef¬ 
fect of radiation feedback can be studied in galaxies that 
evolve in their natural environment, and we can study the 
effect on galaxy evolution, inflows, outflows, and the observ¬ 
able properties of the ISM and GGM. Going to smaller scales 
will allow us to properly resolve optically thick star-forming 
clouds and how they are affected by stellar. 
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APPENDIX A: STELLAR LUMINOSITIES AND 
PHOTON PROPERTIES 


The emission from each stellar particle is calculated on the 
fly for every hne timestep and injected into the host cell, 
adding to the radiation energy density of all photon groups. 
For the specific stellar luminosity (i.e. luminosity per unit 
mass) and photon group properties, we use the spectral 


emission distribution (SED) models of Bruzual & Chariot 


(2003), where we assume a Chabrier (2003) IMF. The de¬ 
pendence of the specific luminosities and radiation group 
properties on the stellar population’s age and metallicity 
are shown in Fig. A1 Each photon group’s properties (i.e. 
average energy and cross section) are updated every 5 coarse 
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Figure Al. Age and metallicity dependence of specific stellar luminosities and photon group attributes derived from the |Bruzual fc| 
[Chariot] ( |2003| > SED model, assuming a jChabrierj ( |200^ IMF. The columns represent the five photon groups with increasing photon 
energy from left to right. Top row: specific luminosity (i.e. luminosity per unit stellar mass), emitted into each photon group from the 
stellar particles. Second row: cumulative specific luminosity per photon group. Third row: average photon energies. Bottom row: 
average cross sections per photoionisation interaction. The emission from each stellar particle is calculated on-the-fiy in each timestep 
by integration of the data shown in the second row, given the mass, age and metallicity of the particle. The properties of the five photon 
groups are updated every five coarse timesteps by a luminosity-weighted average of all existing stellar particles (excluding the stellar 
particles present in the initial conditions). 


timesteps, using luminosity-weighted averages of the exist¬ 
ing stellar particles’ emission in the corresponding bands us¬ 
ing the frequency dependent ionisation cross sections from 
Verner et al. ( 1996| (see also the appendix of R13). This up¬ 


date of the photon groups and the stellar emission is done 
as detailed in the appendix of R13, except that the specific 
luminosity is now in terms of emitted energy whereas it was 
done in terms of photon number count in R13 (which makes 
more sense in pure ionisation calculations). The stellar emis¬ 
sion is thus energy-conserving, whereas it was photon num¬ 
ber conserving in R13. This difference arises because the 
spectral shape of an individual stellar particle is not identi¬ 
cal to the “average” shape which is assumed for our photon 
groups. Due to this difference, one must choose whether the 
emission is accurate in terms of photon count or energy, and 
we have chosen energy. 


APPENDIX B: THE REDUCED FLUX 
APPROXIMATION 

In ^ we describe the reduced flux approximation, whereby 
we assume a full reduced flux of photons, |F| = cE^ when 


calculating the direct radiation pressure on gas from non- 
IR photons. The reason for making this approximation is 
as follows. Radiation is emitted from a stellar particle di¬ 
rectly into the cell which hosts the particle, by increment¬ 
ing the radiation energy density, while the photon flux 
is left unchanged, in accordance with locally isotropic ra¬ 
diation from the stellar population. We use the so-called 
Global Lax Friedrich Riemann (GLF) solver for the advec- 
tion of photons between cells (see R13, R14), which has the 
advantage that radiation is advected isotropically from such 
sources, i.e. the radiation field retains an isotropic shape, 
in the limit of a single source and free-streaming radiation. 
The disadvantage of the GLF solver is that the radiation 
stays somewhat isotropic inside a buffer of a few cell widths 
around such a source, i.e. in this region |F| <C cE. 

We demonstrate this in Fig. |B1| where we show the 
converged results of a 3-d experiment of a single isotropic 
source of radiation in the middle of a box resolved by 64^ 
cells, and assuming free-streaming radiation, i.e. no interac¬ 
tion between the radiation and the medium. The photons are 
injected by incrementing E during each timestep, uniformly 
in the eight cells adjacent to the box center, according to the 
luminosity of the source. We plot, as a function of distance 
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Figure Bl. Converged results from a simple 3-d experiment of 
a single isotropic source of free-streaming radiation in the center 
of a box resolved by 64^ cells. We plot the analytic expectation 
for the r~^ radiation flux (dashed blue), the “angle-integrated” 
radiation flux, cE (solid red), and the magnitude of the radiation 
flux away from the source, |F| (solid green), all in units of Fq, 
the expected analytic flux at a distance of one cell width from 
the source. Against the right axis, in purple, we plot the reduced 
flux, /^, which is the ratio of |F| and cE. The reduced flux curve 
demonstrates a “smoothing length” for the reduced flux of a few 
cell widths around the isotropic radiation source, within which 
the radiation pressure, oc F, is considerably underestimated. 
The curves for cF, |F|, and have been binned by radius, and 
the thickness of the curves reflects the flux range within each 
radial bin. 

r from the source, the analytic expression for the radiation 
flux, i.e. the “angle-integrated” flux cF, and the magnitude 
of the photon bulk flux, |F| (pointing away from the source), 
all in units of Fq = ^^^^2 ? where L is the source luminosity 
and Ax the cell width. It can be seen from the plot that 
cE follows the analytic profile accurately, but within a few 
Ax from the source, |F| ^ cE. In solid purple, against the 
right axis, we plot also the reduced flux of the radiation, 
f'y ~ While in reality, one would have = 1 at any ra¬ 
dius for this simple experiment, this is clearly quite far from 
the truth close to the source, with e.g. f^{r < bAx) < 0.8. 
For the advection of photons, photoionisation and the asso¬ 
ciated heating, this is of no consequence, but the radiation 
pressure is correspondingly underestimated in such a buffer 
of ~ 5 cell widths around the source, which can be consid¬ 
ered a “smoothing length” for the radiation pressure oc F 
(see R14, Sec. 2.3.3, and Eq. 27). Typical Hii regions in our 
simulations are badly resolved, which means that most of 
the ionising radiation is absorbed within 5 cell widths from 
the stellar sources, and hence the direct radiation pressure 
is potentially underestimated. 

We have therefore used the aforementioned reduced flux 
approximation, re-normalising F to cE for the radiation 
pressure force in each cell, i.e. assuming a full reduced flux, 
in the bulk direction of the radiation. As discussed in ^ this 
means we overestimate the radiation pressure in two types of 
locations: i) In cells hosting stellar radiation sources, where 
the radiation is in reality isotropic, but we instead take it 
all to point in the same (average) direction, ii) In-between 
radiation sources, where the radiation pressure from oppos¬ 
ing fields of radiation would in reality cancel out, but again 



Figure Cl. Resolution tests. The plot shows comparisons of the 
stellar mass formed in the g9 galaxy for default (thick curves) 
and low (thin curves) resolution, for the various feedback models, 
as indicated in the legend. 



Figure C2. Resolution tests with constant stellar particle mass 
(600 Mq). The plot shows comparisons of the stellar mass formed 
in the g9 galaxy for default (thick curves) and low (thin curves) 
resolution, for the various feedback models, as indicated in the 
legend. 

we instead take it to point in the average direction, which is 
the direction away from the strongest source. Since we found 
the effect of direct radiation pressure to be negligible, the 
use of the reduced flux approximation is conservative and 
our conclusions are robust. 


APPENDIX C: RESOLUTION TESTS 

As is the case with simulation work in general, it is im¬ 
portant to investigate the dependence of the results on the 
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numerical resolution. For this reason, we compare runs of 
the fiducial g9 galaxy with lower-resolution counterparts, 
where the minimum cell width is two times larger and the 
mass of particles (both formed and in the initial conditions) 
is 8 times largei|^ Otherwise the simulation parameters are 
identical. 

Fig. |C1| shows the effect of resolution on star forma¬ 
tion, with line thickness indicating the resolution (thick for 
high resolution and thin for low resolution), and as usual 
the colours and linestyles represent different feedback mod¬ 
els. We skip here the results from the runs comparing the 
different radiation processes, i.e. including/excluding direct 
and reprocessed radiation pressure, and show only the all- 
inclusive radiation runs, but note that the radiation pressure 
processes have little and seemingly random effects, i.e. the 
dominant radiation effect is heating, as we have established 
in the previous sections. 

Lowering the resolution has the effect of reducing the 
formation of stars, regardless of the feedback process in¬ 
cluded, even without any feedback. This is a natural out¬ 
come of lowering the resolution, since it becomes more diffi¬ 
cult for the gas to collapse to high densities, which in turn 
decreases the star formation rate which scales locally as 
Indeed we find the mean densities typically to be higher in 
the high-resolution runs than in their low-resolution coun¬ 
terparts, by about half a dex. Apart from this systematic 
suppression in star formation rates with resolution, stellar 
radiation feedback reduces the formation of stars by roughly 
a similar fraction: without SN feedback, the addition of ra¬ 
diation feedback reduces the stellar mass formed at 500 Myr 
by 40%. Combined SN and RT feedback results in very 
similar star formation for the two resolutions, indicating nu¬ 
merical convergence. 

In the resolution comparison, it is questionable whether 
m* should be changed in the low resolution simulations. In¬ 
creasing the stellar particle mass by a factor of 8, as we have 
done in Fig. |C1| can boost the feedback, since each parti¬ 
cle then has eight times higher luminosity and SN energy, 
thus contaminating the pure effect of changing the physical 
resolution. For the sake of completeness, we thus also ran 
lower-resolution counterparts to the g9 simulations, exactly 
as just described, but with m* = 600 M©, identical to the 
fiducial simulations, for which the star formation is com¬ 
pared to the higher-resolution case in Fig. |C2| Here the ef¬ 
fect of SN feedback is negligible at low resolution, while pure 
radiation feedback is more efficient than shown in Fig. 


affect our results. For this purpose, we have run the lower- 
resolution equivalent of the g9 galax^Qwith light speeds six 
times lower, two times higher, and six times higher than the 
default value, the last value being recommended by |Rosdahl| 


et al. (2013). 


In Fig. m we plot a comparison of the different light 
speed runs in the form of the total stellar mass formed during 
the 500 Myr run time. The light speed has a negligible effect 
on the star formation. The morphology, outflow rates, and 
density distributions are also nearly identical in the different 
light speed runs. 

We conclude that our results are well converged in terms 
of the employed light speed, and we expect that similar re¬ 
sults would be retrieved with the full light speecQ 

While true in the main simulation runs described in 
this paper, our conclusion on light speed convergence does 
not necessarily hold when the IR optical depth becomes 
very high, as in our “extreme” simulations described in §4.3| 
where we artificially boosted the IR opacity by orders of 
magnitude compared to the more realistic theoretically mo¬ 
tivated value, and found very reduced star formation and 
outflows. When the optical depth becomes very high, the ef¬ 
fective propagation speed of radiation scales inversely with 
the local optical depth, i.e. radiation waves travel at a speed 
c/r, where r is the optical depth across some relevant length 
scale, such as an optically thick cloud (see e.g. sections 2.4 
and 3.5 in R14). With our reduced speed of light, radiation 
waves travel at a speed c/r, and if c is orders of magnitude 
smaller than the real light speed, as in this paper, the speed 
of light can become a severe issue in very optically thick 
gas, with radiation waves potentially travelling at a speed 
slower than the gas itself. Since this becomes most severe 
with the highest optic al depths, we ran the most extreme 
experiment from 14.3 (^ir = lO^cm^g”^) with c increased 
by a factor two and decreased by a factor three from the 
default value, i.e. c = 10“^ c and c = 1.67 x 10“^ c, re¬ 
spectively. We found that the average star formation rates 
are unaffected, but that they fluctuate on a longer timescale 
with decreasing light speed, and that outflow rates decrease 
very substantially with increasing light speed (reinforcing 
our conclusion that radiation does not produce outflows). It 
thus appears that with very large optical depth, the reduced 
speed of light does become an issue for outflows, but not for 
star formation rates. The main conclusions of this paper are 
not affected though, since optical depths in our simulations 
are never very high, save for the extreme “what if” scenario 
described in 84.31 


APPENDIX D: REDUCED LIGHT SPEED 
CONVERGENCE TESTS 


To prevent a prohibitively small timestep in our RHD 
scheme, we use a default reduced light speed of c = c/200 
in our simulation runs. This is in fact a six times lower light 
speed than recommended in the analysis of reduced light 
speeds in ISM simulations in Rosdahl et al. (2013), so it 
is important to verify that the chosen light speed does not 


® Correspondingly, the initial condition particles are 8 times 
fewer. 


^ eight times fewer/more massive particles and twice the min¬ 
imum cell width compared to the default resolution, just as in 
Appendix [C] 

® Apart from the problem of the computational cost of a run with 
the full light speed, such a simulation would also likely suffer from 
hydrodynamical diffusion with the current setup, due to a large 
number of very small timesteps. For a full light speed to work 
with our explicit RT solver, we would need to sub-cycle the RT 
within the HD step. 
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Figure Dl. Light speed convergence tests. The plot shows com¬ 
parisons of the stellar mass formed in the g 9 lower resolution 
galaxy (see Sec. for the default light speed of c/200 (solid 
green), along with identical runs with the light speed changed by 
factors of 6, 2, and 1/3, as indicated in the legend. 
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